/* OpenTeddy 3d sketch modeler based on the paper by T. Igarashi et al.
*      "Teddy: A Sketching Interface for 3D Freeform Design"
*
* Copyright (c) 2004 Alex Ni
*/


#include "otconfig.h"

#include "GTS/gts.h"

#include <math.h>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <time.h>
#include <string>

extern "C"
{
#include "triangle.h"
#include "trackball.h"
}

// C++
#include <map>
#include <list>
#include <deque>
#include <vector>

// Our includes
#include "main.h"
#include "pub_infl.h"

#include "../glfw/win32/glfw.h"

#include "otstu.h"
#include "stuff.h"


double g_scale = 1.0;
double g_scale_actual = 1.0;
int win_w;
int win_h;
double proj_w;
double proj_h;

ot_public_inflate* g_opi;


//hack
double ot_depth_range()
{
	return 20.0;
}

void ot_set_projection()
{
	if (win_w > win_h)
	{
		proj_h = 1.0;
		proj_w = ((double)win_w) / ((double)win_h);
	}
	else
	{
		proj_w = 1.0;
		proj_h = ((double)win_h) / ((double)win_w);
	}
	glOrtho( -proj_w, proj_w, -proj_h, proj_h, -(ot_depth_range() / 2.0), ot_depth_range() / 2.0);
	//gluPerspective(45.0, proj_w/proj_h, 1.0, 100.0);
}
// setup opengl viewing matrices
void setup_viewing_matrices()
{
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	ot_set_projection();

	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
}

bool manifold_warning = false;
bool closed_warning = false;

void ot_set_geometry_warning(bool mw, bool cw)
{
	manifold_warning = mw;
	closed_warning = cw;
}

void ot_calculate_geometry_warning(GtsSurface* s)
{
	ot_set_geometry_warning(
		!gts_surface_is_manifold(s),
		!gts_surface_is_closed(s)
	);
}

static char g_string_to_display[256] = {0};

void ot_set_display_string(char* s, bool log = true)
{
	if (strlen(s) > 255) return ;
	strcpy(g_string_to_display, s);
	g_string_to_display[(strlen(s))] = '\0';

	if (log)
		printf("Display: %s \n", s);

	g_string_to_display[0] = 0;
}
void ot_set_display_string_modal(char* s)
{
	ot_set_display_string(s, false);
}

void ot_clear_display_string()
{
	g_string_to_display[0] = '\0';
}

void ot_display_status()
{
	char to_display[512];
	sprintf(to_display, "%s%s%s",
			manifold_warning ? "M " : "",
			closed_warning ? "C " : "",
			g_string_to_display);
	printf(to_display);
}


double ot_loop_perimeter(OT2loop* points, bool closed)
{
	double total_length = 0.0;
	for (int i = 0;i < (closed ? points->size() : (points->size() - 1));i++)
	{
		total_length += (*points)[i].dist( (*points).wrap(i + 1) );
	}

	return total_length;
}

#define OT_MIN(a,b) ( ((a)<(b)) ? (a) : (b) )

// place points at regular intervals
void ot_loop_regulate(OT2loop* points, bool closed, int dnp)
{
	OT2loop pnew(*points);
	double l = ot_loop_perimeter(points, closed);

	// desired length of each segment
	if (!closed) dnp--;
	double dlen = l / dnp;

	// re-push_back desired points into points
	points->clear();


	double dleft = dlen;  // distance left until next point
	double dedge = pnew[0].dist(pnew[1]);  // distance left on current edge
	double dedgel = dedge; // total current edge length
	points->push_back( pnew[0] );

	int pnewi = 0; // beginning point index of current edge
	dnp--; // since we did the beginning point, don't do as many other points

	for (int i = 0; i < dnp; i++)
	{
		dleft = dlen;
		while (dleft > 0.0)
		{ // if we still have distance to go before next pt
			if (dleft > dedge)
			{ // edge is shorter than distance to next pt
				dleft -= dedge;

				pnewi++; // go to next edge
				dedgel = dedge = pnew.wrap(pnewi).dist(pnew.wrap(pnewi + 1));
			}
			else
			{ // edge still has some distance on it
				dedge -= dleft;

				dleft = 0.0;
			}
		}

		// plot a new point
		points->push_back( pnew.wrap(pnewi).interpolate(pnew.wrap(pnewi + 1), 1.0 - dedge / dedgel) );
	}
	if (!closed) points->push_back( pnew[pnew.size() - 1] );
}

void ot_loop_simplify(OT2loop* points, bool closed, double granularscale = 1.0)
{
	OT_ASSERT(points);

	if (points->size() < 10) return ;  // sort of an arbitrary number, may want to change

	double l = ot_loop_perimeter(points, closed);

	// desired number of points
	int dnp = l / (granularscale * ot_get_pref_double("loop_simp_dist", 0.04)) - 1; // at this point, arbitrary

if (dnp < 10){ dnp = OT_MIN(10, points->size()); }

	ot_loop_regulate(points, closed, dnp);
}

bool ot_segment_intersects_polygon(OT2vec v1, OT2vec v2, OT2loop* points, bool do_exclude_index = false, size_t exclude_index = 0)
{
	bool intersection = false;

	GtsVertex* p1, *p2, *pnew1, *pnew2;
	GtsSegment* e_new, *e_old;
	pnew1 = gts_vertex_new (gts_vertex_class (),
			v2[0], v2[1], 0.0 );
	pnew2 = gts_vertex_new (gts_vertex_class (),
			v1[0], v1[1], 0.0 );

	p1 = gts_vertex_new (gts_vertex_class (),
			0.0, 0.0, 0.0 );
	p2 = gts_vertex_new (gts_vertex_class (),
			0.0, 0.0, 0.0 );
	e_new = gts_segment_new( gts_segment_class(), pnew1, pnew2 );
	e_old = gts_segment_new( gts_segment_class(), p1, p2 );
	for (size_t i = 0;i < points->size();i++)
	{
		if (do_exclude_index && (i == exclude_index)) continue;
		p1->p.x = (*points)[i][0];
		p1->p.y = (*points)[i][1];
		p2->p.x = (*points).wrap(i + 1)[0];
		p2->p.y = (*points).wrap(i + 1)[1];
		if (GTS_IN == gts_segments_are_intersecting(e_new, e_old))
		{
			intersection = true;
		}
	}

	// this should take care of the points
	gts_object_destroy (GTS_OBJECT (e_new));
	gts_object_destroy (GTS_OBJECT (e_old));

	return intersection;
}

#define MAX_PTS 1024
#define PT(i,j)		(*points)[i][j]

void addpt(double x, double y, OT2loop* points, bool intersect_okay = false)
{
	size_t ptct = points->size();
	if (ptct >= MAX_PTS) return ;

	if ( !intersect_okay && (ptct > 3) && ot_segment_intersects_polygon((*points)[ptct - 1], OT2vec(x, y), points))
	{
		points->clear();
	}
	else
	{
		points->push_back( OT2vec(x, y) );
	}
}

bool ot_polygon_intersecting(OT2loop* points)
{
	for (int i = 0;i < points->size();i++)
	{
		if (ot_segment_intersects_polygon( (*points)[i], (*points).wrap(i + 1), points, true, i ))
		{
			return true;
		}
	}
	return false;
}

double polygon_area(OT2loop* points)
{
	double total_area = 0.0;
	size_t ptct = points->size();
	for (int i = 0;i < ptct;i++)
	{
		// cross product
		total_area += (*points)[i].cross( (*points).wrap(i + 1) );
	}
	return total_area;
}



static double mvmatrix[16], projmatrix[16];
static GLint viewport[4];

void update_matrices()
{
	glGetDoublev(GL_MODELVIEW_MATRIX, mvmatrix);
	glGetDoublev(GL_PROJECTION_MATRIX, projmatrix);
	glGetIntegerv(GL_VIEWPORT, viewport);
}

void mproj(double x, double y, double z, double* rx, double *ry, double *rz)
{
	gluProject(x, y, z, mvmatrix, projmatrix, viewport,
			rx, ry, rz);
}
double screen_to_input(double sc, double proj, double win)
{
	return (2*proj*sc / win) - proj;
}
void mproj_input(double x, double y, double z, double* rx, double *ry, double *rz)
{
	gluProject(x, y, z, mvmatrix, projmatrix, viewport,
			rx, ry, rz);
	*rx = screen_to_input(*rx, proj_w, win_w);
	*ry = screen_to_input(*ry, proj_h, win_h);
}

void getmouse(int* x = NULL, int* y = NULL, double* dx = NULL, double* dy = NULL)
{
	int ix, iy;
	glfwGetMousePos(&ix, &iy);
	iy = win_h - iy;

	if (x) *x = ix;
	if (y) *y = iy;
	if (dx) *dx = screen_to_input(ix, proj_w, win_w);
	if (dy) *dy = screen_to_input(iy, proj_h, win_h);
}


void geteyevector(double *x, double *y, double *z)
{
	// note, gl projection default is looking down the negative z axis.
	// hence, project v1=(0,0,1) and v2=(0,0,0) and v2-v1 gives eye direction
	double v1[3], v2[3];
	mproj(0.0, 0.0, 1.0, &v1[0], &v1[1], &v1[2]);
	mproj(0.0, 0.0, 0.0, &v2[0], &v2[1], &v2[2] );
	for (int i = 0;i < 3;i++){ v2[i] -= v1[i]; }

	gts_vector_normalize(v2);
	*x = v2[0]; *y = v2[1]; *y = v2[2];
}

void munproj(double x, double y, double d, double* rx, double *ry, double *rz, bool no_proj = false)
{
	if (no_proj)
	{
		//GLdouble old_w = projmatrix[10];
		//projmatrix[10] = 0.5;

		gluUnProject(x, y, d,
				mvmatrix, projmatrix, viewport,
				rx, ry, rz);

		//projmatrix[10] = old_w;
	}
	else
	{
		gluUnProject(x, y, d,
				mvmatrix, projmatrix, viewport,
				rx, ry, rz);
	}
}

double input_to_screen(double in, double proj, double win)
{
	return ((in + proj) / (2*proj))*win;
}

void munproj_input(double x, double y, double d, double* rx, double *ry, double *rz, bool no_proj = false)
{
	munproj( input_to_screen(x, proj_w, win_w), input_to_screen(y, proj_h, win_h), d,
			rx, ry, rz, no_proj );
}

bool faces_front(GtsFace* f)
{
	GtsVertex *v[3];
	gts_triangle_vertices(&f->triangle,
			&v[0], &v[1], &v[2]);

	OT3vec win[3];
	for (int i = 0;i < 3;i++)
	{
		mproj(v[i]->p.x, v[i]->p.y, v[i]->p.z,
				&win[i][0], &win[i][1], &win[i][2]);
		win[i][2] = 0.0;
	}
	// calculate cross product
	win[1] -= win[0];
	win[2] -= win[0];
	OT3vec xprod;
	xprod.set_cross(win[1], win[2]);

	return xprod[2] > 0.0;
}

// in global coordinates, what is the x dimension at depth d?
static double ot_measure_screen_at_depth(double d)
{
	OT3vec c1, c2;
	munproj_input(0.0, 1.0, d, &c1[0], &c1[1], &c1[2]);
	munproj_input(1.0, 1.0, d, &c2[0], &c2[1], &c2[2]);
	c1 -= c2;
	return c1.length();
}



// change coord from -proj_w, proj_w,  proj_h, -proj_h to
static void ot_gts_vertex_to_global(GtsVertex* v, double* d)
{
	double sd;
	sd = ot_measure_screen_at_depth(*d);

	OT3vec newv;

	munproj_input(v->p.x, v->p.y, *d + (1.0 / g_scale_actual) * -(1.0 / ot_depth_range())*(sd * v->p.z), &newv[0], &newv[1], &newv[2], true);
	v->p.x = newv[0]; v->p.y = newv[1]; v->p.z = newv[2];
}
static void ot_gts_vertices_to_global(GtsSurface* s, double d)
{
	gts_surface_foreach_vertex(s, (GtsFunc) ot_gts_vertex_to_global, &d);
}

// flip all faces of a gts surface
//////////////////
static void flip_face(GtsFace* f, GtsSurface *s)
{
	gts_triangle_revert(&f->triangle);
}

static void flip_faces(GtsSurface* gtss)
{
	gts_surface_foreach_face(gtss, (GtsFunc) flip_face, gtss);
}
///////////////////

typedef std::vector<GtsSegment*> seglist_t;
typedef std::vector<GtsVertex*> vertlist_t;

GtsVertex* nextsegvert(GtsSegment* s, GtsVertex* cur)
{
	if (cur == s->v1) return s->v2;
	if (cur == s->v2) return s->v1;
	OT_ASSERT(0);
	return NULL;
}

// this part of GTS should be made public
static gint triangle_orientation (GtsTriangle * t, GtsEdge * e)
{
	GtsSegment * s = GTS_SEGMENT (t->e1 == e ? t->e2
			:
			t->e2 == e ? t->e3
			:
			t->e1);
	GtsVertex * v2 = GTS_SEGMENT (e)->v2;

	if (s->v1 == v2 || s->v2 == v2)
		return 1;
	return -1;
}



void makevertlist(seglist_t* segl, vertlist_t* vl, GtsSurface* sur)
{
	seglist_t::iterator sli = segl->begin();

	// safe to do this because all of our segments are edges
	GtsVertex *v = NULL;
	{
		GSList singleedge; singleedge.data = ((GtsEdge*)(*sli)); singleedge.next = NULL;
		GSList* faces = gts_faces_from_edges(&singleedge, sur);
		OT_ASSERT(faces && faces->data);
		GtsFace* f = (GtsFace*)faces->data;
		if (triangle_orientation((GtsTriangle*)f, (GtsEdge*)(*sli)))
			v = (*sli)->v1;
		else
			v = (*sli)->v2;
	}

	OT_ASSERT(v != NULL);


	GtsVertex *i = v;
	do
	{
		vl->push_back(i);
		i = nextsegvert( (*sli), i );
		sli++;
	}
	while ( sli != segl->end() );

	//vl now has all points of segl
}

GtsEdge* ot_gts_find_edge_or_new(GtsSurface* s, GtsVertex* v1, GtsVertex* v2)
{
	OT_ASSERT(v1 != NULL); OT_ASSERT(v2 != NULL);

	if (v1 == v2) return NULL;

	GtsEdge* ee = NULL;

	// create constant 2 item list
	GSList vlist2; vlist2.data = v2; vlist2.next = NULL;
	GSList vlist; vlist.data = v1; vlist.next = &vlist2;

	GSList* edges = gts_edges_from_vertices(&vlist2, s);
	GSList* l = edges;
	while (l)
	{
		GtsEdge* e = ((GtsEdge*)l->data);
		if (	(e->segment.v1 == v1 && e->segment.v2 == v2) ||
				(e->segment.v1 == v2 && e->segment.v2 == v1) )
		{
			ee = e;
		}
		l = l->next;
	}
	g_slist_free(edges);

	if (ee == NULL) ee = gts_edge_new( s->edge_class, v1, v2 );

	OT_ASSERT(ee != NULL);
	return ee;
}

// v0
// |  \
// e2   e0
// |       \
// v2 - e1 - v1
//
GtsFace* ot_gts_add_tri_from_existing_verts(GtsSurface* s, GtsVertex** v )
{
	GtsFace* f;
	GtsEdge* e[3];
	for (int i = 0;i < 3;i++)
	{
		e[i] = ot_gts_find_edge_or_new(s, v[i], v[(i + 1) % 3]);

		OT_ASSERT(e[i]);
	}
	gts_surface_add_face (s, f = gts_face_new (s->face_class,
			e[0], e[1], e[2] ) );

	for (int i = 0;i < 3;i++)
	{
		OT_ASSERT(e[i]->triangles != NULL);
	}
	return f;
}

// v0 - e0 - v1
// |  \       |
// e3   dia   e1
// |       \  |
// v2 - e2 - v3
//
// note, dia is not checked in all cases- may need to add that in the future
// if the mechanism is needed (e.g. suppose there already exists an edge
// connecting v2 and v1, or possibly both v0-v3 and v1-v2, then we cannot
// do anything

bool ot_gts_add_quad_from_existing_verts(GtsSurface* s, GtsVertex** v )
{
	// null quad cases
	if ( (v[0] == v[3]) || (v[2] == v[1]) ) return false;

	if ( (v[0] == v[1]) && (v[2] == v[3]) ) return false;
	if ( (v[2] == v[0]) && (v[1] == v[3]) ) return false;

	// triangle cases
	if ( (v[0] == v[1]) || (v[1] == v[3]) )
	{
		GtsVertex* vt[3];
		vt[0] = v[0];
		vt[1] = v[3];
		vt[2] = v[2];
		ot_gts_add_tri_from_existing_verts(s, vt );
		return true;
	}
	if ( (v[0] == v[2]) || (v[2] == v[3]) )
	{
		GtsVertex* vt[3];
		vt[0] = v[0];
		vt[1] = v[1];
		vt[2] = v[3];
		ot_gts_add_tri_from_existing_verts(s, vt );
		return true;
	}


	GtsEdge* e[4];
	GtsEdge* dia;

	int choice[] = {0, 1, 3, 2};

	for (int i = 0;i < 4;i++)
	{
		e[i] = ot_gts_find_edge_or_new(s, v[choice[i]], v[choice[(i + 1) % 4]]);
		OT_ASSERT(e[i]);
	}

	// create triangles
	dia = ot_gts_find_edge_or_new(s, v[0], v[3] );
	gts_surface_add_face (s, gts_face_new (s->face_class,
			e[0], e[1], dia ) );

	gts_surface_add_face (s, gts_face_new (s->face_class,
			e[2], e[3], dia ) );

	return true;
}

// NOTE- INPUT WORKING CORRECTLY DEPENDS IF POINTS ARE IN CLOCKWISE OR COUNTERCLOCKWISE ORDER
// MAKE SURE YOU KNOW WHICH IS WHICH

// extrude the points in 'points' towards viewing direction.
// if nearpts is not NULL, add all near points to it (used to extract closest loop)
// if tiponly is true, only make the tip part (NOT YET IMPLEMENTED, right now whole thing is always made)
GtsSurface* generate_slice(OT2loop* points, std::list<GtsVertex*> *nearpts, bool tiponly)
{
	GtsSurface* slice_s = gts_surface_new (gts_surface_class (),
			GTS_FACE_CLASS (gts_face_class ()),
			GTS_EDGE_CLASS (gts_edge_class ()),
			GTS_VERTEX_CLASS (ot_gts_vertex_class ()));

	GtsVertex* slicepts[4];
	GtsEdge* sliceedges[4];
	GtsFace* sliceface[2];
	GtsSegment* sliceseg[2];
	GtsEdge* first_sliceedge = NULL;
	GtsVertex* first_slicepts[2];
	first_slicepts[0] = first_slicepts[1] = NULL;

	GtsVertex* tips[2];
	GtsEdge* tipedges[2][2];
	GtsEdge* first_tipedge[2];
	first_tipedge[0] = first_tipedge[1] = NULL;

	tips[0] = tips[1] = NULL;
	tipedges[0][0] = tipedges[0][1] = tipedges[1][0] = tipedges[1][1] = NULL;

	int ptct = points->size();
	for (int ptn = 0;ptn < ptct;ptn++)
	{
		// try to cut trhough
		for (int i = 0;i < 4;i++)
		{
			if ( (i >= 2) && (ptn == (ptct - 1) ) )
			{
				slicepts[i] = first_slicepts[i - 2];
				// wrap around
				continue;
			}

			bool makenew = true;
			if (i < 2)
			{ // try to use the last ones if possible
				if (ptn > 0)
				{
					slicepts[i] = slicepts[2 + i];
					makenew = false;
				}
			}
			if (makenew)
			{
				slicepts[i] = gts_vertex_new( slice_s->vertex_class, 0.0, 0.0, 0.0);
				if ( nearpts && (!tiponly) && ((i % 2) == 0) ) nearpts->push_back(slicepts[i]);
				munproj_input( (*points).wrap(ptn + i / 2)[0], (*points).wrap(ptn + i / 2)[1],
						tiponly ? ((i % 2) ? 100.0 : 1.0) : ((i % 2) ? 1.0 : 0.0),  //1.0 for far, 0.0 for near
								&slicepts[i]->p.x, &slicepts[i]->p.y, &slicepts[i]->p.z);
				// save first two
				if ( !first_slicepts[(i % 2)] ) first_slicepts[i % 2] = slicepts[i];
			}
		}
		// create tip points of cylinder (making it 'pointy', like 1 cone on each side)
		// in tiponly mode we only want the near cone.
		if (!tips[0])
		{
			for (int i = 0;i < 2;i++)
			{
				tips[i] = gts_vertex_new( slice_s->vertex_class, 0.0, 0.0, 0.0);
				if (tiponly && nearpts && (i == 0)) nearpts->push_back(tips[i]);
				munproj( win_w / 2, win_h / 2,
						tiponly ? ((i % 2) ? 100.0 : -100.0) : ((i % 2) ? 1.0 : 0.0),   // -100 very near
								&tips[i]->p.x, &tips[i]->p.y, &tips[i]->p.z);
			}
		}


		GtsEdge* diagonal;
		int choice[] = {0, 1, 3, 2};
		for (int i = 0;i < 4;i++)
		{
			if ( (i == 2) && (ptn == (ptct - 1) ) )
			{
				// wrap around
				sliceedges[i] = first_sliceedge;
			}
			else
			{
				if ( (ptn > 0) && (i == 0) ) sliceedges[i] = sliceedges[2];
				else
					sliceedges[i] = gts_edge_new( slice_s->edge_class, slicepts[choice[i]], slicepts[choice[(i + 1) % 4]]);
			}

			if (!first_sliceedge) first_sliceedge = sliceedges[i];
		}

		// slice triangles
		diagonal = gts_edge_new( slice_s->edge_class, slicepts[0], slicepts[3] );
		gts_surface_add_face (slice_s, sliceface[0] = gts_face_new (slice_s->face_class,
				sliceedges[0], sliceedges[1], diagonal ) );

		gts_surface_add_face (slice_s, sliceface[1] = gts_face_new (slice_s->face_class,
				sliceedges[2], sliceedges[3], diagonal ) );
		// tip triangles (cone)

		// for far and near
		for (int i = 0;i < 2;i++)
		{
			for (int j = 0;j < 2;j++)
			{
				if ( (j == 1) && (ptn == (ptct - 1) ) )
				{
					// wrap around
					tipedges[i][j] = first_tipedge[i];
				}
				else
				{
					if ( (ptn > 0) && (j == 0) ) tipedges[i][j] = tipedges[i][1];
					else
						tipedges[i][j] = gts_edge_new( slice_s->edge_class, tips[i], slicepts[i ? (j * 2 + 1) : (j * 2)] );
				}
				if (!first_tipedge[i]) first_tipedge[i] = tipedges[i][j];
			}
		}

		for (int i = 0;i < 2;i++)
		{
			for (int j = 0;j < 2;j++)
			{
				OT_ASSERT(tipedges[i][j] != NULL);
			}
		}

		gts_surface_add_face (slice_s, gts_face_new(slice_s->face_class, tipedges[0][0], sliceedges[3], tipedges[0][1] ) ); // far
		gts_surface_add_face (slice_s, gts_face_new(slice_s->face_class, tipedges[1][0], tipedges[1][1], sliceedges[1] ) ); // near
		// end tip


		// slice rays
		sliceseg[0] = gts_segment_new( gts_segment_class(),
				slicepts[2], slicepts[0] );
		sliceseg[1] = gts_segment_new( gts_segment_class(),
				slicepts[3], slicepts[1] );

	}

	// face right way
	if (polygon_area(points) < 0.0 )
	{
		flip_faces(slice_s);
	}

	return slice_s;
}

#define MARK_KEEP ( (gpointer) 1)
#define MARK_THROW_AWAY ( (gpointer) 2)

static void mark_gts_object (GtsObject * o, gpointer data)
{
	o->reserved = data;
}
static void mark_gts_object_and (GtsObject * o, gpointer data)
{
	if (o->reserved) o->reserved = data;
	else o->reserved = NULL;
}

static guint remove_marked(GtsFace* f, GtsSurface *s)
{
	return (GTS_OBJECT(f)->reserved == MARK_THROW_AWAY);
}

GtsSurface* ot_gts_boolean_union(GtsSurface* s1, GtsSurface* s2)
{
	bool success = true;

	GNode* s1tree = gts_bb_tree_surface(s1);
	GNode* s2tree = gts_bb_tree_surface(s2);
	GtsSurfaceInter* si = NULL;
	GtsSurface* gtss = NULL;

	if (!s1tree){ success = false; goto ot_gts_boolean_union_done; }
	if (!s2tree){ success = false; goto ot_gts_boolean_union_done; }

	// new resulting surface
	gtss = gts_surface_new (gts_surface_class (),
			GTS_FACE_CLASS (gts_face_class ()),
			GTS_EDGE_CLASS (gts_edge_class ()),
			GTS_VERTEX_CLASS (ot_gts_vertex_class ()));

	si = gts_surface_inter_new (gts_surface_inter_class (),
			s1, s2, s1tree, s2tree, FALSE, FALSE);
	if (!si){ success = false; goto ot_gts_boolean_union_done; }
	gboolean isclosed;
	gts_surface_inter_check(si, &isclosed);
	if (!isclosed){ success = false; goto ot_gts_boolean_union_done; }

	/* This is the correct command for boolean union */
	gts_surface_inter_boolean (si, gtss, GTS_2_OUT_1);
	gts_surface_inter_boolean (si, gtss, GTS_1_OUT_2);

ot_gts_boolean_union_done:
	// destroy all temporary objects
	if (s1tree) gts_bb_tree_destroy(s1tree, TRUE);
	if (s2tree) gts_bb_tree_destroy(s2tree, TRUE);

	if (si) gts_object_destroy (GTS_OBJECT (si));

	if (success)
		return gtss;
	else
	{
		gts_object_destroy (GTS_OBJECT (gtss) );
		return NULL;
	}
}


#define OT_OBJ_RESERVED(v) (GTS_OBJECT(v)->reserved)


// transform vertex by matrix
static void ot_gts_vertex_matmult(GtsVertex* gtsv, double* m)
{
	OT3vec v; v[0] = gtsv->p.x; v[1] = gtsv->p.y; v[2] = gtsv->p.z;
	v.matmult(m);
	gtsv->p.x = v[0]; gtsv->p.y = v[1]; gtsv->p.z = v[2];
}
static void ot_gts_vertices_matmult(GtsSurface* s, double* m)
{
	gts_surface_foreach_vertex(s, (GtsFunc) ot_gts_vertex_matmult, m);
}

// transform vertex by matrix (once only)
static void ot_gts_vertex_matmult_norepeat(GtsVertex* gtsv, double* m)
{
	if (OT_OBJ_RESERVED(gtsv)) return ;
	OT3vec v; v[0] = gtsv->p.x; v[1] = gtsv->p.y; v[2] = gtsv->p.z;
	v.matmult(m);
	gtsv->p.x = v[0]; gtsv->p.y = v[1]; gtsv->p.z = v[2];
	OT_OBJ_RESERVED(gtsv) = (gpointer)true;
}
static void ot_gts_vertices_matmult_norepeat(GtsSurface* s, double* m)
{
	gts_surface_foreach_vertex(s, (GtsFunc) ot_gts_vertex_matmult_norepeat, m);
}

////////////////////////////////
//  drawing stuff
////////////////////////////////

// drawing

enum render_mode_t { OT_RENDER_WIREFRAME = 0, OT_RENDER_SOLID };

render_mode_t g_render_mode = OT_RENDER_WIREFRAME;

bool g_show_wireframe = true;


OT3vec g_eyevec;
// main gts drawing interface
inline static void draw_face_simple(GtsFace*f, OT2vec* loc)
{
	GtsVertex *v[3];
	gts_triangle_vertices(&f->triangle,
			&v[0], &v[1], &v[2]);
	glVertex3dv( &v[0]->p.x );
	glVertex3dv( &v[1]->p.x );
	glVertex3dv( &v[2]->p.x );
}

inline static void draw_face(GtsFace* f, OT2vec* loc)
{
	GtsFace *of = GTS_FACE(f);

	GtsVertex *v[3];
	gts_triangle_vertices(&f->triangle,
			&v[0], &v[1], &v[2]);

	for (int i = 0;i < 3;i++)
	{
		if (g_render_mode != OT_RENDER_WIREFRAME)
		{
			OTGtsVertex *ov = OT_GTS_VERTEX(v[i]);
			double dot = -((ov->n).dot(g_eyevec)); //dot=abs(dot);
			if (dot < 0.0) dot = 0.0;

			glColor3f( dot, dot, dot);
		}
		//glColor3f(1.0, 0.0, 0.0);
		glVertex3dv( &v[i]->p.x );
		//glVertex3d( v[i]->p.x, v[i]->p.y, v[i]->p.z);
	}
	// the above might fuck up if Gts implementation changes, but it's faster for now
	//		glVertex3d( v[i]->p.x, v[i]->p.y, v[i]->p.z);
	//	}
}

////////////////////////////
static bool g_auto_select;

OT2loop g_debugloop; //TEST

typedef std::map<GtsVertex*, double> OTvpdhash;

class OTblob
{
public:
	bool invis;
	bool dirty_facesel;
	bool bSelected;
	GtsSurface* gtss;
	GtsSurfaceInter* si;
	vertlist_t traceverts;
	OTvpdhash vertsel;
	OTvpdhash vertlabel;

	// ordering on the faces
	std::vector<GtsFace*> indf;
	// ordering on the vertices and keyframe delta transformed vertices
	std::vector<GtsVertex*> indv;


	OTblob()
	{
		invis = false;

		dirty_facesel = true;

		bSelected = g_auto_select;
		gtss = NULL;
		si = NULL;
	}
	~OTblob()
	{
		clear();
	}
	void clear_si()
	{
		if (si) gts_object_destroy(GTS_OBJECT(si));
		si = NULL;
	}

	void clear_gtss()
	{
		if (gtss) gts_object_destroy (GTS_OBJECT (gtss));
		gtss = NULL;
	}
	void matmult(double* mat)
	{
		OT_ASSERT(mat);

		if (gtss) ot_gts_vertices_matmult_norepeat(gtss, mat);
		if (si)
		{
			if (si->s1) ot_gts_vertices_matmult_norepeat(si->s1, mat);
			if (si->s2) ot_gts_vertices_matmult_norepeat(si->s2, mat);
		}
		if (gtss) gts_surface_foreach_vertex(gtss, (GtsFunc) gts_object_reset_reserved, NULL);
		if (si)
		{
			if (si->s1) gts_surface_foreach_vertex(si->s1, (GtsFunc) gts_object_reset_reserved, NULL);
			if (si->s2) gts_surface_foreach_vertex(si->s2, (GtsFunc) gts_object_reset_reserved, NULL);
		}
	}

	// create an ordering on the faces
	void index_faces();

	GtsFace* pick_closest_face(int x, int y, bool unsel, bool laser = false)
	{
		index_faces();
		if (indf.size() == 0) return NULL;

		GLint viewport[4];
		glGetIntegerv(GL_VIEWPORT, viewport);

		glMatrixMode(GL_PROJECTION);
		glPushMatrix();
		glLoadIdentity();
		gluPickMatrix(x, y, 1.0, 1.0, viewport); //(win_w+win_h)/40, (win_w+win_h)/40, viewport);
		ot_set_projection();

		if (laser)
			glDisable(GL_CULL_FACE);
		else
			glEnable(GL_CULL_FACE);

		GLuint *selectBuf = new GLuint[indf.size() * 2];
		glSelectBuffer(indf.size()*2, selectBuf);
		glRenderMode(GL_SELECT);
		glInitNames();
		glPushName(0);

		for (int i = 0;i < indf.size();i++)
		{
			glLoadName(i);
			glBegin(GL_TRIANGLES);
			draw_face(indf[i], NULL);
			glEnd();
		}
		glPopMatrix();
		glFlush();


		GLint hits = glRenderMode(GL_RENDER);

		GLuint *ptr = selectBuf;
		for (int i = 0;i < hits;i++)
		{
			GLuint names = *ptr++;
			// skip z
			ptr++; ptr++;
			for (int j = 0;j < names;j++)
			{
				int newind = *ptr++;
				if (newind >= indf.size() || newind < 0) continue;

				GtsFace* of = GTS_FACE(indf[newind]);

				GtsVertex *v[3];
				gts_triangle_vertices(&indf[newind]->triangle,
						&v[0], &v[1], &v[2]);
				for (int i = 0;i < 3;i++)
					unsel ? unsetsel(v[i]) : select(v[i], 1.0);

			}
		}

		delete[] selectBuf;
		glEnable(GL_CULL_FACE);

		return NULL;
	}

	void move_vertex(GtsVertex *v, OT3vec loc)
	{
		v->p.x = loc[0]; v->p.y = loc[1]; v->p.z = loc[2];
	}

	void select(GtsVertex* v, double d)
	{
		vertsel[v] = d;
		dirty_facesel = true;
	}

	double getsel(GtsVertex* v)
	{
		OTvpdhash::iterator i = vertsel.find(v);
		if (i == vertsel.end()) return 0.0;
		return (*i).second;
	}

	void unsetsel(GtsVertex* v)
	{
		OTvpdhash::iterator i = vertsel.find(v);
		if (i != vertsel.end()) vertsel.erase(i);
	}

	void clearsel()
	{
		vertsel.clear();
		dirty_facesel = true;
	}

	void selsel();
	void comp_normals();

	void clear()
	{
		clear_si(); clear_gtss();
		traceverts.clear();
		clearsel();
	}
	bool join(OTblob* b)
	{
		GtsSurface* gtss_union =
			ot_gts_boolean_union(gtss, b->gtss);
		if (!gtss_union) return false;
		clear();
		gtss = gtss_union;
		return true;
	}
};


// create an ordering on the faces of the surface
static void ot_gts_add_face_to_index(GtsFace * f, OTblob* b)
{
	b->indf.push_back(f);
}

void OTblob::index_faces()
{
	// check if presumably we have already done it
	//if( gts_surface_face_number(gtss) == indf.size() ) return;

	indf.clear();
	gts_surface_foreach_face (gtss, (GtsFunc) ot_gts_add_face_to_index, this);
}


static void ot_gts_calc_normal(GtsVertex *v, OTblob* b)
{
	static OTGtsVertex *ov;ov = OT_GTS_VERTEX(v);

	/* loop subd... won't work right since neighbors are not ordered
	// calc tangent vectors
	OT3vec t1(0.0,0.0,0.0), t2(0.0, 0.0, 0.0);

	GSList* allneigh = gts_vertex_neighbors( v, NULL, b->gtss );

	{
		GSList* lis = allneigh;
		OTdouble len = g_slist_length(allneigh);
		OTdouble i = 0.0;
		while(lis){
			OTdouble weight = cos( (2*M_PI)*(i/len) );
			t1[0] += weight*GTS_VERTEX(lis->data)->p.x;
			t1[1] += weight*GTS_VERTEX(lis->data)->p.y;
			t1[2] += weight*GTS_VERTEX(lis->data)->p.z;

			weight = sin( (2*M_PI)*(i/len) );
			t2[0] += weight*GTS_VERTEX(lis->data)->p.x;
			t2[1] += weight*GTS_VERTEX(lis->data)->p.y;
			t2[2] += weight*GTS_VERTEX(lis->data)->p.z;

			i += 1.0;

			lis = lis->next;
		}
	}
	g_slist_free(allneigh);

	ov->n.set_cross(t1, t2);
	ov->n.normalize();
	*/

	static GSList* allface;allface = gts_vertex_faces( v, b->gtss, NULL );

	ov->n = OT3vec(0.0, 0.0, 0.0);
	{
		static GSList* lis;lis = allface;
		//OTdouble len = g_slist_length(allneigh);
		//OTdouble i = 0.0;
		while (lis)
		{
			static OT3vec cn;
			gts_triangle_normal(GTS_TRIANGLE(lis->data), &cn[0], &cn[1], &cn[2]);
			ov->n += cn;

			//i += 1.0;

			lis = lis->next;
		}
	}
	g_slist_free(allface);

	ov->n.normalize();

	/*
		if(gts_vertex_mean_curvature_normal(v, b->gtss, &ov->n[0])){
			ov->n.normalize();
			// probably should have it point out the right way
		}
	//	v->n = OT3vec(x,y,z);
	*/
}

void OTblob::comp_normals()
{
	gts_surface_foreach_vertex (gtss, (GtsFunc) ot_gts_calc_normal, this);
}

static void ot_gts_add_surround_sel(GtsVertex * v, OTblob* b)
{

	double dist_to_s1 = 1e50; // infinity would be nice

	for (int i = 0;i < 2;i++)
	{
		GSList* j = v->segments;
		while (j)
		{
			GtsVertex *adjv;
			GtsSegment* s1 = (GtsSegment*)j->data;

			adjv = NULL;
			if (s1->v1 == v) adjv = s1->v2;
			else adjv = s1->v1;
			GtsVector del;
			gts_vector_init(del, &v->p, &adjv->p);
			double dell = OT3vec(del[0], del[1], del[2]).length();

			if (i == 0)
			{
				// part of first surface
				if ( GTS_OBJECT(adjv)->reserved == (gpointer)1 )
				{
					if (dell < dist_to_s1) dist_to_s1 = dell;
				}
			}
			else
			{
				// not part of first surface
				if ( GTS_OBJECT(adjv)->reserved == (gpointer)2 )
				{
					b->select(adjv, dist_to_s1 / (dist_to_s1 + dell));
				}
			}
			j = j->next;
		}
	}
}
static void ot_gts_add_sel(GtsVertex * v, OTblob* b)
{
	if ( GTS_OBJECT(v)->reserved == (gpointer)2 )
		b->select(v, 1.0);
}
static void ot_gts_force_sel(GtsVertex * v, OTblob* b)
{
	b->select(v, 1.0);
}

void OTblob::selsel()
{
	if (!si)
	{ // select all/unselect all
		if (vertsel.size() > 0)
		{
			clearsel();
		}
		else
		{
			gts_surface_foreach_vertex (gtss, (GtsFunc) ot_gts_force_sel, this);
		}
		return ;
	}
	clearsel();
	gts_surface_foreach_vertex (si->s1, (GtsFunc) gts_object_reset_reserved, NULL);
	gts_surface_foreach_vertex (si->s2, (GtsFunc) gts_object_reset_reserved, NULL);
	gts_surface_foreach_vertex (gtss, (GtsFunc) mark_gts_object, (gpointer)2);

	// part of surface that definitely has
	GtsSurface* sel_gtss = gts_surface_new (gts_surface_class (),
			GTS_FACE_CLASS (gts_face_class ()),
			GTS_EDGE_CLASS (gts_edge_class ()),
			GTS_VERTEX_CLASS (ot_gts_vertex_class ()));
	gts_surface_inter_boolean (si, sel_gtss, GTS_1_IN_2);
	gts_surface_foreach_vertex (sel_gtss, (GtsFunc) ot_gts_add_sel, this);
	gts_surface_foreach_vertex (sel_gtss, (GtsFunc) mark_gts_object, (gpointer)1);


	//gts_surface_foreach_vertex (sel_gtss, (GtsFunc) mark_gts_object, (gpointer)1);

	// follow intersection loop
	for (vertlist_t::iterator i = traceverts.begin();
			i != traceverts.end(); ++i)
	{
		ot_gts_add_surround_sel((*i), this);
	}

	gts_surface_foreach_vertex (sel_gtss, (GtsFunc) gts_object_reset_reserved, NULL);

	gts_object_destroy (GTS_OBJECT (sel_gtss) );
	gts_surface_foreach_vertex (gtss, (GtsFunc) gts_object_reset_reserved, NULL);
	// this is no longer valid
	clear_si();
	traceverts.clear();
}


void ot_gts_new_origin(OTblob* b)
{
	double mat[16];
	glGetDoublev(GL_MODELVIEW_MATRIX, mat);
	/*	for(int i=0;i<16;i++){
			mat[i] = 0.5 * ((i/4) == (i%4));
		}
	*/
	b->matmult(mat);
}


enum perform_action_t { OT_PERFORM_PUNCHOUT = 0, OT_PERFORM_CUT, OT_PERFORM_INV_CUT, OT_PERFORM_PERFORATE};
/* NOTE TO SELF- WILL NEED TO KEEP ORIGINAL GTSS AROUND TO UNDO ACTION */
bool perform_action(perform_action_t which, OTblob* b)
{
	OT_ASSERT(b);
	if (!b->si) return false;

	GtsSurface* oldgtss = b->gtss;

	// new resulting surface
	GtsSurface* gtss = gts_surface_new (gts_surface_class (),
			GTS_FACE_CLASS (gts_face_class ()),
			GTS_EDGE_CLASS (gts_edge_class ()),
			GTS_VERTEX_CLASS (ot_gts_vertex_class ()));

	if (which == OT_PERFORM_PUNCHOUT)
	{

		// add new tesselated to inter_gtss
		//	gts_surface_inter_boolean (si, inter_gtss, GTS_1_IN_2);   // TEMP
		//	gts_surface_inter_boolean(si, inter_gtss, GTS_2_OUT_1);
		gts_surface_inter_boolean (b->si, gtss, GTS_1_OUT_2);
		// inter_Gtss now has the right thing in it.
	}
	else
	{  // CUT
		// do a cut out
		//gts_surface_inter_boolean (si, inter_gtss, GTS_1_OUT_2);
		//gts_surface_inter_boolean (si, inter_gtss, GTS_1_IN_2);
		//gts_surface_inter_boolean (si, inter_gtss, GTS_2_OUT_1);  /*cylinder w/ surface cut out*/
		//gts_surface_inter_boolean (si, inter_gtss, GTS_2_IN_1);   /* cylinder and surface intersection */

		if (which == OT_PERFORM_CUT)
		{
			/* This will give the sliced area cut out */
			gts_surface_inter_boolean (b->si, gtss, GTS_2_IN_1);   /* cylinder and surface intersection */
			flip_faces(gtss);
			gts_surface_inter_boolean (b->si, gtss, GTS_1_OUT_2);
		}
		else if (which == OT_PERFORM_INV_CUT)
		{ // inverse cut
			gts_surface_inter_boolean (b->si, gtss, GTS_2_IN_1);   /* cylinder and surface intersection */
			gts_surface_inter_boolean (b->si, gtss, GTS_1_IN_2);
		}
		else if (which == OT_PERFORM_PERFORATE)
		{ // add edges where we cut thru
			gts_surface_inter_boolean (b->si, gtss, GTS_1_OUT_2);
			gts_surface_inter_boolean (b->si, gtss, GTS_1_IN_2);
		}
	}

	int a = gts_surface_is_closed(gtss);
	int be = gts_surface_is_manifold(gtss);
	printf("Perform action %d... closed:  %d,   mainfold:   %d\n", which, a, be);

	// don't need these any more
	b->clear_si();
	gts_object_destroy (GTS_OBJECT (oldgtss));

	// set new object to gtss
	b->gtss = gtss;

	return true;
}

static bool ot_realize_loop_on_surface(OT2loop* points, vertlist_t* traceverts, OTblob* b, bool do_intersect)
{
	OT_ASSERT(b);
	OT_ASSERT(traceverts);
	OT_ASSERT(points);
	if (points->size() < 3 ) return false;

	bool status = true;

	seglist_t tracesegs;

	b->si = NULL;

	GNode* modeltree = NULL;
	GNode* slicetree = NULL;

	std::list<GtsVertex*> nearpts;
	GtsSurface* slice_s = generate_slice(points, &nearpts, true /*tiponly*/);
	GtsSurface* gtss = b->gtss;
	//gtss = slice_s; return;// TEEEMP
	{
		int a = gts_surface_is_closed(slice_s);
		int b = gts_surface_is_manifold(slice_s);
		printf("Slice... closed:  %d,   mainfold:   %d\n", a, b);
		OT_ASSERT(a);
		OT_ASSERT(b);
	}

	slicetree = gts_bb_tree_surface(slice_s);
	modeltree = gts_bb_tree_surface(gtss);

	b->si = gts_surface_inter_new (gts_surface_inter_class (),
			gtss, slice_s, modeltree, slicetree, FALSE, FALSE);
	// now, we want to go and get the closest (to the camera) loop, since we only want to use that one

	// make sure the cut is valid (should be, provided both slice and gtss are closed manifolds)
	gboolean isclosed;
	gts_surface_inter_check(b->si, &isclosed);
	if (!isclosed || do_intersect)
	{
		if (!do_intersect)
		{
			if (!isclosed)
			{
				status = false;
			}
		}
		else
		{
			if (b->si->edges != NULL)
			{
				status = true;
			}
			else
			{
				status = false;
			}
		}
		if (status == false) b->clear_si();
	}
	else
	{

		// adjust the 3d intersection so that the only intersection
		// we actually use is the closest one (one that has vertices adjacent to nearpts)

#define MARK_VERTEX_NEAR ((gpointer) 1)

		// mark all vertices next to nearpts in slice part of b->si (b->si->s2)
		for (std::list<GtsVertex*>::iterator i = nearpts.begin(); i != nearpts.end(); i++)
		{
			GSList *neighbors = gts_vertex_neighbors( (*i), NULL, b->si->s2 );
			g_slist_foreach(neighbors, (GFunc)mark_gts_object, MARK_VERTEX_NEAR);
			g_slist_free(neighbors);
		}

#define EDGE_NOT_FOLLOWED ((gpointer) 1)
#define EDGE_FOLLOWED ((gpointer) 0)

		typedef std::list<GtsSegment*> edgelist_t;
		typedef std::list<edgelist_t*> edgelistlist_t;
		edgelistlist_t eloops;
		edgelist_t* nearloop = NULL;  // eventually will be set to the closest

		GSList* i = b->si->edges;
		g_slist_foreach (b->si->edges, (GFunc) mark_gts_object, EDGE_NOT_FOLLOWED);
		while (i)
		{
			GtsSegment* seg = (GtsSegment*)i->data;

			if ( (GTS_OBJECT (seg)->reserved) == EDGE_NOT_FOLLOWED )
			{ // not yet taken, follow this edge.
				edgelist_t* curloop = new edgelist_t();
				eloops.push_back(curloop);
				curloop->push_back(seg);			// push back the very first segment in this loop.

				while (GTS_OBJECT (seg)->reserved == EDGE_NOT_FOLLOWED )
				{
					GTS_OBJECT (seg)->reserved = EDGE_FOLLOWED;

					// of all the segments of each vertex, see if we can follow.
					GSList* j;
					for (int i = 0;i < 2;i++)
					{
						GtsVertex* v = (i ? seg->v1 : seg->v2);
						if ( !nearloop && GTS_OBJECT(v)->reserved == MARK_VERTEX_NEAR )
						{
							// we want this loop, since it's adjacent to the near slice pts
							nearloop = curloop;
						}

						j = v->segments;
						while (j)
						{
							GtsSegment* s1 = (GtsSegment*)j->data;
							if ( (s1 != seg) && (GTS_OBJECT(s1)->reserved == EDGE_NOT_FOLLOWED) )
							{
								// follow this new edge
								curloop->push_back(s1);
								seg = s1;
								goto foundedge;
							}

							j = j->next;
						}
					}
foundedge:
					;
				}

			}
			i = i->next;
		}
		g_slist_foreach (b->si->edges, (GFunc) gts_object_reset_reserved, NULL);


		// we want to keep only the closest loop, start with marking to throw away all triangles
		//#define MARK_KEEP ( (gpointer) 1)
		//#define MARK_THROW_AWAY ( (gpointer) 1)

		gts_surface_foreach_face (b->si->s1, (GtsFunc) mark_gts_object, MARK_THROW_AWAY);


		// follow all loops
		int x, y;
		x = g_slist_length(b->si->edges);
		int counter = 0;
		for (edgelistlist_t::iterator i = eloops.begin(); i != eloops.end(); i++)
		{
			if (nearloop == (*i)) tracesegs.clear();

			for (edgelist_t::iterator j = (*i)->begin(); j != (*i)->end(); j++)
			{
				if (nearloop == (*i)) // for the near loop, we want to keep
				{
					tracesegs.push_back((*j));

					/* MARK ALL TRIANGLES OF EACH VERTEX OF (*j) to "KEEP" */

					GSList* allfaces = gts_vertex_faces( (*j)->v1, b->si->s1, NULL );
					allfaces = gts_vertex_faces( (*j)->v2, b->si->s1, allfaces );

					GSList* lis = allfaces;
					while (lis)
					{
						GTS_OBJECT((GtsFace*)lis->data)->reserved = MARK_KEEP;
						lis = lis->next;
					}
					g_slist_free(allfaces);
				} else
				{
					// this could be more efficient...
					// remove all edges not  in near loop
					b->si->edges = g_slist_remove(b->si->edges, (*j));

					/* MARK ALL TRIANGLES OF EACH VERTEX OF (*j) to "THROW AWAY" IF THEY ARE NOT "KEEP" */
				}
			}

			counter++;
			// delete this loop
			delete (*i);
		}
		/* THROW AWAY ALL TRIANGLES IN b->si->s1 THAT ARE NOT "KEEP" */
		{
			gts_surface_foreach_face_remove (b->si->s1, (GtsFunc) remove_marked, b->si->s1);
		}

		gts_surface_foreach_face (b->si->s1, (GtsFunc) gts_object_reset_reserved, NULL);

		// unmark all vertices next to nearpts in slice part of b->si (b->si->s2)
		for (std::list<GtsVertex*>::iterator i = nearpts.begin(); i != nearpts.end(); i++)
		{
			GSList *neighbors = gts_vertex_neighbors( (*i), NULL, b->si->s2 );
			g_slist_foreach(neighbors, (GFunc)gts_object_reset_reserved, NULL);
			g_slist_free(neighbors);
		}

		//gtss = si->s1; return; //TEMP_
		/* temp, check number of edges has changed */

		y = g_slist_length(b->si->edges);

		printf("%d %d\n", x, y);
		/* temp check number of edges has changed */

		makevertlist(&tracesegs, traceverts, b->si->s1);
	} // if(!closed)


	//
	// destroy all temporary objects
	gts_bb_tree_destroy(modeltree, TRUE);
	gts_bb_tree_destroy(slicetree, TRUE);

	// destroy original surfaces (since we now have the edges from them we want)

	gts_object_destroy (GTS_OBJECT (slice_s));

	b->gtss = gtss;
	return status;
}



////////////////////////////////
//  drawing stuff
////////////////////////////////

GtsFace *g_closest_face;
GLfloat g_zerodepth;
GLfloat g_olddepth;

static void draw_face_only(GtsFace* f, OTblob* b)
{
	glBegin(GL_TRIANGLES);
	draw_face(f, NULL);
	glEnd();
}

int countit = 0;

int olcount;
static void outline_face2(GtsFace* f, GtsSurface *s)
{
	if (olcount++ != countit) return ;
	glBegin(GL_LINE_LOOP);
	draw_face_simple(f, NULL);
	glEnd();
}

double g_lastalp;
static void ot_gts_draw_sel_vert(GtsVertex* v, OTblob* b)
{
	double alp;
	if ( (alp = b->getsel(v)) <= 0.0) return ;

	// can buy a little speed from not sending a color down the graphics pipeline every time
	if (alp != g_lastalp)
	{
		glColor4d(0.0, 1.0, 0.0, alp);
		g_lastalp = alp;
	}

	//	glColor4d(alp, 1.0-alp, 0.0, 1.0); // temp

	glBegin(GL_POINTS);
	glVertex3dv( &v->p.x );
	// the above might fuck up if Gts implementation changes, but it's faster for now
	//		glVertex3d( v->p.x, v->p.y, v->p.z);
	//	}

	glEnd();
}

static void outline_face(GtsFace* f, GtsSurface *s)
{
	glBegin(GL_LINE_LOOP);
	draw_face_simple(f, NULL);
	glEnd();
}
////////////////////////////////
// end more drawing stuff
////////////////////////////////




inline void normalize_2d(double* x, double* y)
{
	double len = ot_2d_dist(*x, *y, 0, 0);
	*x /= len; *y /= len;
}

inline double dotprod_2d(double x0, double y0, double x1, double y1)
{
	return x0*x1 + y0*y1;
}

// find the tangent at each point, return the sum of dot product with x & y
double tangentdot(OT2loop* points, GLuint p1, GLuint p2)
{
	int ptct = points->size();
	// make chord
	double lx = input_to_screen(PT(p1, 0), proj_w, win_w);
	double rx = input_to_screen(PT(p2, 0), proj_w, win_w);
	double ly = input_to_screen(PT(p1, 1), proj_h, win_h);
	double ry = input_to_screen(PT(p2, 1), proj_h, win_h);

	double cseclen = ot_2d_dist(lx, ly, rx, ry );
	double x = (lx - rx) / cseclen; double y = (ly - ry) / cseclen;

	// done

	OT2vec tgt[2];
	GLint choices[4]; // which points to choose

	// first tangent
	choices[0] = p1 - 1;
	choices[1] = p1 + 1;
	// second tangent
	choices[2] = p2 + 1;
	choices[3] = p2 - 1;

	for (int i = 0;i < 4;i++)
	{
		if (choices[i] < 0) choices[i] = 0;
		if (choices[i] >= ptct) choices[i] = ptct - 1;
	}

	// calc tangents, normalize
	for (int i = 0;i < 2;i++)
	{
		for (int j = 0;j < 2;j++)
		{
			tgt[i][j] = input_to_screen(PT(choices[i * 2 + 1], j), (j == 0) ? proj_w : proj_h, (j == 0) ? win_w : win_h )
					- input_to_screen(PT(choices[i * 2], j), (j == 0) ? proj_w : proj_h, (j == 0) ? win_w : win_h ) ;
		}

		if ( (tgt[i][0] == 0.0) && (tgt[i][1] == 0.0) )
		{
			printf("Warning: zero-length tangent\n");
			return 0.0; // can't really do anything else at this point
		}
		normalize_2d(&tgt[i][0], &tgt[i][1]);
	}

	double result = 0.0;
	double results[2];
	for (int i = 0;i < 2;i++)
	{
		//result += abs(M_PI/2 - acos(abs(dotprod_2d(tgt[i][0], tgt[i][1], x, y))));
		result += abs(dotprod_2d(tgt[i][0], tgt[i][1], x, y));
		results[i] = abs(dotprod_2d(tgt[i][0], tgt[i][1], x, y));
		results[i] = abs((M_PI / 2) - acos(results[i]));
	}


	// average func
	//return results[0]*0.5 + results[1] * 0.5;

	// min func
	if (results[0] < results[1]) return results[0];
	return results[1];

	return result;
}

#define EXTRUDE_CLOSE_ENOUGH 0.17

static bool extrude_trace(OT2loop* points, vertlist_t* traceverts, OTblob* b)
{
	OT_ASSERT(b);

	if (points->size() <= 0) return false;
	if (!traceverts) return false;
	// special case, if we only have two points, we assume the user wants a flat surface
	if (points->size() == 2)
	{
		ot_loop_regulate(points, false, 3);
	}

	int ptct = points->size();

	// find centroid of loop
	double vavg[3];
	{
		GLuint cnt = 0;
		for (int i = 0;i < 3;i++) vavg[i] = 0.0;
		for (vertlist_t::iterator i = traceverts->begin();
				i != traceverts->end(); ++i)
		{
			vavg[0] += (*i)->p.x; vavg[1] += (*i)->p.y; vavg[2] += (*i)->p.z;
			cnt++;
		}
		for (int i = 0;i < 3;i++) vavg[i] /= cnt;
	}
	// find what depth coordinate the centroid would have, we use that one
	double davg;
	{
		double win[2];
		mproj(vavg[0], vavg[1], vavg[2], &win[0], &win[1], &davg);
	}

	double p0[2], pn[2];
	p0[0] = input_to_screen(PT(0, 0), proj_w, win_w); p0[1] = input_to_screen(PT(0, 1), proj_h, win_h);
	pn[0] = input_to_screen(PT(ptct - 1, 0), proj_w, win_w); pn[1] = input_to_screen(PT(ptct - 1, 1), proj_h, win_h);

	double pd[2], pc[2]; // delta, center
	for (int i = 0;i < 2;i++){ pd[i] = pn[i] - p0[i]; pc[i] = (pn[i] + p0[i]) / 2.0; }

	// normalize
	double plen = ot_2d_dist(pd[0], pd[1], 0, 0);
	pd[0] /= plen; pd[1] /= plen;


	// ok, here's what's easier to do
	// for each point on the loop
	//   project the point pL (find where it is currently in screen space)
	//   take the dot product (pL-p0) * pd, this is the new 'y' coordinate
	//     (the y coordinate is where the loop point is placed on the current slice of the extrusion)

	double* ycoord = new double[traceverts->size()];
	double* depth = new double[traceverts->size()];
	double depthavg = 0.0;
	{
		double maxy;
		double miny;

		{
			int count = 0;
			for (vertlist_t::iterator i = traceverts->begin();
					i != traceverts->end(); ++i)
			{
				// project tracevert
				double win[3];
				mproj((*i)->p.x, (*i)->p.y, (*i)->p.z, &win[0], &win[1], &win[2]);
				// dot 2d loc with pd
				ycoord[count] = (win[0] - p0[0]) * pd[0] + (win[1] - p0[1]) * pd[1];
				if (count == 0){ maxy = miny = ycoord[count]; }
				else
				{
					if (ycoord[count] > maxy) maxy = ycoord[count];
					if (ycoord[count] < miny) miny = ycoord[count];
				}

				depthavg += (depth[count] = win[2]);
				count++;
			}
			depthavg /= count;
		}


		// normalize ycoords  to lie between [0,1]
		for (int i = 0;i < traceverts->size();i++)
		{
			ycoord[i] -= miny;
			ycoord[i] = ycoord[i] / (maxy - miny);

			// center depth on depthavg
			depth[i] -= depthavg;
		}

	}

	// build up new verts, stitch together
	// start with our original traceverts, stitch up to first and last points on stroke
	GtsVertex** lastpts = new GtsVertex * [traceverts->size()];
	GtsVertex** curpts = new GtsVertex * [traceverts->size()];

	{
		int count = 0;
		for (vertlist_t::iterator i = traceverts->begin();
				i != traceverts->end(); ++i)
		{
			curpts[count] = NULL;
			lastpts[count++] = (*i);
		}
	}

	int new_lowpts = 1;
	int new_highpts = 1;  // high or low is advanced, make new points

	size_t countend = traceverts->size();

	size_t left = 0; size_t right = ptct - 1;
	bool done = (right - left) < 2;  // special case if we have only one point (in which case, make a 'cone')
	if (done)
	{	// make the traceverts be the final ring
		left = -1;
		for (int i = 0;i < countend;i++)
		{
			curpts[i] = lastpts[i];
		}
	}
	while ( !done )
	{
		// get new cross section length
		double lx = input_to_screen(PT(left, 0), proj_w, win_w);
		double rx = input_to_screen(PT(right, 0), proj_w, win_w);
		double ly = input_to_screen(PT(left, 1), proj_h, win_h);
		double ry = input_to_screen(PT(right, 1), proj_h, win_h);

		double cseclen = ot_2d_dist(lx, ly, rx, ry );
		double csecdx = (lx - rx) / cseclen; double csecdy = (ly - ry) / cseclen;


		{
			size_t count = 0;
			for (count = 0; count < countend; count++)
			{
				double win[3];
				win[0] = (1.0 - ycoord[count]) * lx + (ycoord[count]) * rx;
				win[1] = (1.0 - ycoord[count]) * ly + (ycoord[count]) * ry;
				win[2] = (depth[count] * (cseclen / plen)) + depthavg;
				// project this win point in, add it
				GtsVertex *v = NULL;
				if (!new_highpts && (ycoord[count] == 1.0) )
				{ // Use previous 'high' point if any
					v = lastpts[count];
				}
				if (!new_lowpts && (ycoord[count] == 0.0) )
				{ // Use previous 'low' point if any
					v = lastpts[count];
				}
				if (!v)
				{	// else, make a new point
					v = gts_vertex_new( b->gtss->vertex_class, 0.0, 0.0, 0.0);
					munproj(win[0], win[1], win[2],
							&v->p.x, &v->p.y, &v->p.z);
					traceverts->push_back(v);
				}
				curpts[count] = v;
			}
			// stitch curpts to lastpts
			for (count = 0; count < countend; count++)
			{
				GtsVertex* vs[4];
				vs[0] = curpts[count]; vs[1] = curpts[ (count + 1) % countend ];
				vs[2] = lastpts[count]; vs[3] = lastpts[ (count + 1) % countend ];

				ot_gts_add_quad_from_existing_verts(b->gtss, vs);
			}

			// move cur points to last points
			for (count = 0; count < countend; count++)
			{
				lastpts[count] = curpts[count];
			}
			// move cur edges to last edges
		}

		// the fact that there isn't a good way to just say, here's some stuff, pick the min in C++ is sad

		// choose best way to advance points

		int choices[3][2] = { {1, 0}, {0, -1}, {1, -1} };
		{
			int bestchoice = 0; double min;
			double fchoices[3];
			double cur;

			for (int i = 0;i < 3;i++)
			{
				if ( (i == 2) && ((right - left) < 4) ) continue;
				cur = abs(tangentdot(points, left + choices[i][0], right + choices[i][1]));
				if (i == 0) min = cur;
			else if (cur <= min){ bestchoice = i; min = cur; }
				fchoices[i] = cur;
			}
			printf("%d %f %f %f\n", bestchoice, fchoices[0], fchoices[1], fchoices[2]);
			// if close enough, we prefer to advance both
			if ( ( (right - left) >= 4 ) && ( abs(min - fchoices[2]) < EXTRUDE_CLOSE_ENOUGH ) ) bestchoice = 2;



			if ( (right - left) > 2)
			{ // if it's the last time, don't increment left and right
				new_lowpts = choices[bestchoice][0];
				left += choices[bestchoice][0];

				new_highpts = choices[bestchoice][1];
				right += choices[bestchoice][1];
			}
			else
			{
				done = true;
			}
		}
	}

	// now, connect the last ring to the final point.
	GtsVertex *v = gts_vertex_new( b->gtss->vertex_class, 0.0, 0.0, 0.0);
	munproj_input(PT(left + 1, 0), PT(left + 1, 1), depthavg,
			&v->p.x, &v->p.y, &v->p.z);
	traceverts->push_back(v);

	for (int count = 0; count < countend; count++)
	{
		GtsVertex* vs[3];
		vs[1] = curpts[count]; vs[0] = curpts[ (count + 1) % countend ];
		vs[2] = v;

		ot_gts_add_tri_from_existing_verts(b->gtss, vs);
	}

	delete[] lastpts;
	delete[] curpts;
	// for each point on the loop
	//   now we have a location on the slice of the extrusion, scale the depth of pL based on length
	//  ... ok, now we have a 2d point and a depth.  unproject into object space
	//done with one slice of the extrusion!
	delete[] ycoord;
	delete[] depth;

	int a = gts_surface_is_closed(b->gtss);
	int be = gts_surface_is_manifold(b->gtss);
	printf("Extrusion... closed:  %d,   mainfold:   %d\n", a, be);
	ot_calculate_geometry_warning(b->gtss);
	//OT_ASSERT(a);
	//OT_ASSERT(b);

	return true;
}


float g_quat[4], g_spin_quat[4];

static void clear_g_spin_quat()
{
	trackball(g_spin_quat, 0.0, 0.0, 0.0, 0.0);
}
static void clear_g_quat()
{
	trackball(g_quat, 0.0, 0.0, 0.0, 0.0);
}

struct init_g_quats
{
	init_g_quats()
	{
		clear_g_spin_quat();
		clear_g_quat();
	}
}
init_g_quats_singleton;

OT3vec g_cursor, g_newcursor, g_defaultcursor;
bool g_cursor_vis = false;

double rotx, roty;
OT2loop g_points, g_over_points;

void ot_clear_all_modeltrans()
{
	clear_g_spin_quat();
	clear_g_quat();
	rotx = roty = 0.0;
	g_scale = 1.0;
}

typedef std::vector<OTblob*> OTblobvec;
OTblobvec g_allblobs;
OTblob* g_newselsurf = NULL;


void ot_selblob(OTblobvec* ball, OTblobvec* bvsel)
{
	bvsel->clear();
	for (OTblobvec::iterator i = ball->begin(); i != ball->end(); i++)
	{
		OT_ASSERT( (*i) );
		if ( (*i)->bSelected )
		{
			bvsel->push_back( (*i) );
		}
	}
	return ;
}

static double g_global_transform[16];
static double g_global_transform_inv[16];

struct init_g_global_transform
{
	init_g_global_transform()
	{
		for (int i = 0;i < 16;i++){ g_global_transform_inv[i] = g_global_transform[i] = ( (i % 4) == (i / 4) ) ? 1 : 0; }
	}
}
g_global_transform_singleton;

static void ot_modelview_to_global_transform()
{
	glGetDoublev(GL_MODELVIEW_MATRIX, g_global_transform);
	ot_matrix_invert(g_global_transform_inv, g_global_transform);
}

static void ot_add_to_global_transformd(double* mat)
{
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glMultMatrixd(mat);
	glMultMatrixd(g_global_transform);

	ot_modelview_to_global_transform();
}

static void ot_add_to_global_transformf(GLfloat* mat)
{
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glMultMatrixf(mat);
	glMultMatrixd(g_global_transform);

	ot_modelview_to_global_transform();
}

static void ot_view_setup_matrices(bool clear)
{

	if (clear) return ;

	//setglupvector(sin(countit/20.0), cos(countit/20.0) );

	//glRotatef(rotx, 1.0, 0.0, 0.0);
	//glRotatef(roty, 0.0, 1.0, 0.0);
	{
		GLfloat m[4][4];

		build_rotmatrix(m, g_spin_quat);
		glMultMatrixf(&m[0][0]);
	}
	glScalef(1.0 / g_scale, 1.0 / g_scale, 1.0 / g_scale);
}

int draw_2d = 0;


// show a color wheel?
bool g_colorburst = false;
bool g_colorsample = false;
OT3vec g_curcolor(1.0, 1.0, 1.0);

// last place clicked/moved
int eventx = 0;
int eventy = 0;


// if true select face at cursor
int g_pickface = 0;


static void tickdraw()
{
	OTblobvec bsel;
	ot_selblob(&g_allblobs, &bsel);

	glEnable(GL_TEXTURE_2D);
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
	glShadeModel(GL_SMOOTH);

	glDisable(GL_CULL_FACE);
	glDepthFunc(GL_LEQUAL);
	//glEnable(GL_DEPTH_TEST);

	//glEnable(GL_NORMALIZE);


	// draw here
	//if(curlevel) curlevel->draw();
	glClearColor(212.0f / 255.0f, 208.0f / 255.0f, 200.0f / 255.0f, 0.0f);
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

	glDisable(GL_TEXTURE_2D);

	glEnable(GL_CULL_FACE);

	if (draw_2d && g_opi)
	{
		g_opi->draw_2d_debug();
	}

	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	glMultMatrixd(g_global_transform);
	update_matrices();


	int x, y;
	x = eventx; y = eventy;

	if (x < 0) x = 0;
	if (y < 0) y = 0;
	if (y >= win_h) y = win_h - 1;
	if (x >= win_w) x = win_w - 1;

	if (g_allblobs.size() > 0 && !g_colorburst)
	{



		/////////////////////////////// blobs

		int cnt = 0;

		GLfloat olddepth;
		GLfloat newdepth;

		glFlush();
		glReadPixels(x, y, 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &olddepth);
		g_zerodepth = olddepth;
		g_olddepth = g_zerodepth;
		GLuint ndi;
		glReadPixels(x, y, 1, 1, GL_DEPTH_COMPONENT, GL_UNSIGNED_INT, &ndi);
		//printf("Olddepth1: %lf, %u\n", olddepth, ndi);
		olddepth = -0.001f;
		glRasterPos2d(x, y);
		glDrawPixels(1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &olddepth);
		glReadPixels(x, y, 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &olddepth);
		//printf("Olddepth2: %lf\n", olddepth);
		g_newselsurf = NULL;

		munproj(x, y, 0.5, &g_newcursor[0], &g_newcursor[1], &g_newcursor[2]);
		munproj(x, y, 0.5, &g_defaultcursor[0], &g_defaultcursor[1], &g_defaultcursor[2]);
		for (OTblobvec::iterator i = g_allblobs.begin(); i != g_allblobs.end(); i++)
		{
			OT_ASSERT( (*i) );
			if ( (*i)->invis ) continue;

			GtsSurface* gtss = (*i)->gtss;

			if (g_render_mode != OT_RENDER_WIREFRAME) (*i)->comp_normals();

			glMatrixMode(GL_MODELVIEW);
			glLoadIdentity();
			ot_view_setup_matrices( !(*i)->bSelected && (bsel.size() > 0) );
			glMultMatrixd(g_global_transform);

			update_matrices();
			OT3vec centervec;
			munproj(0.0, 0.0, 0.0, &centervec[0], &centervec[1], &centervec[2]);
			munproj(0.0, 0.0, 1.0, &g_eyevec[0], &g_eyevec[1], &g_eyevec[2]);
			g_eyevec -= centervec;
			g_eyevec.normalize();

			// line drawing off
			glDepthMask(GL_TRUE);

			glEnable(GL_DEPTH_TEST);
			if (g_render_mode == OT_RENDER_WIREFRAME)
			{
				ot_set_rgb(( !(*i)->bSelected ) ? ot_get_pref_rgb("bodycolor", OT3vec(1.0, 1.0, 1.0)) :
						ot_get_pref_rgb("selcolor", OT3vec(0.7, 0.7, 1.0)), 1.0);
			}

			glEnable(GL_POLYGON_OFFSET_FILL);
			glPolygonOffset(1.0, 8);


			OT2vec loc(x, y);
			// try geometry intersect -alex
			//ot_click_hit(gtss, input_to_screen(x, proj_w, win_w), input_to_screen(y, proj_h, win_h) );
			// actually this blows, it's too slow


			if (g_pickface && (*i)->bSelected )
			{
				(*i)->pick_closest_face(x, y, g_pickface == 2);
				g_pickface = 0;
			}

			glMatrixMode(GL_MODELVIEW);

			{
				OT2vec loc(x, y);
				if (g_render_mode == OT_RENDER_WIREFRAME)
				{
					glDisable(GL_TEXTURE_2D);
					//glBegin(GL_TRIANGLES);
					gts_surface_foreach_face(gtss, (GtsFunc) draw_face_only, (*i));
					//glEnd();
				}
				else
				{
					glEnable(GL_TEXTURE_2D);
					glTexEnvf(GL_TEXTURE_ENV, GL_TEXTURE_ENV_MODE, GL_MODULATE);
					glShadeModel(GL_SMOOTH);
					gts_surface_foreach_face(gtss, (GtsFunc) draw_face_only, (*i));
				}
			}


			glColor4f(0.0, 0.0, 0.0, 1.00);

			glDisable(GL_TEXTURE_2D);
			// line drawing on,  take off z=writing
			glDepthMask(GL_FALSE);
			glDisable(GL_POLYGON_OFFSET_FILL);
			olcount = 0;
			if (g_show_wireframe)
			{
				gts_surface_foreach_face(gtss, (GtsFunc) outline_face, gtss);
			}

			// DRAW SELECTION POINTS
			// point drawing on
			glDisable(GL_TEXTURE_2D);
			glColor4f(0.0, 1.0, 0.0, 1.0);
			g_lastalp = 1.0;

			glPointSize(3.0);
			glBegin(GL_POINTS);
			gts_surface_foreach_vertex(gtss, (GtsFunc) ot_gts_draw_sel_vert, (*i));
			glEnd();

			// check for new cursor selection
			GLuint ndi;
			glReadPixels(x, y, 1, 1, GL_DEPTH_COMPONENT, GL_UNSIGNED_INT, &ndi);
			glReadPixels(x, y, 1, 1, GL_DEPTH_COMPONENT, GL_FLOAT, &newdepth);

			if (newdepth != olddepth)
			{
				g_newselsurf = (*i);
				olddepth = newdepth;
				munproj((double)x, (double)y, ( (double)newdepth ) / ( (double)g_zerodepth ), &g_newcursor[0], &g_newcursor[1], &g_newcursor[2]);
				//printf("Cursor: %d %d %lf / %lf %lf %lf / %u\n", x, y, (double)newdepth, g_newcursor[0], g_newcursor[1], g_newcursor[2], ndi );
			}

			//glColor4f(0.0, 1.0, 0.0, 1.00);
			//gts_surface_foreach_face(gtss, (GtsFunc) outline_face2, gtss);
			cnt++;

			glPopMatrix();
		}

		glDepthMask(GL_FALSE);
		glDisable(GL_TEXTURE_2D);

		for (OTblobvec::iterator i = g_allblobs.begin(); i != g_allblobs.end(); i++)
		{
			OTblob* b = *i;
			if ( (*i)->invis ) continue;


			/////////////////////////////////// 2sel
			if (b->vertlabel.size() > 1)
			{
				glBegin(GL_LINES);
				glColor3f(1.0, 1.0, 0.0);
				OTvpdhash::iterator i = b->vertlabel.begin();
				OT3vec a( (*i).first->p.x, (*i).first->p.y, (*i).first->p.z );
				i++;
				OT3vec b( (*i).first->p.x, (*i).first->p.y, (*i).first->p.z );
				b -= a;
				OT3vec c = a;
				b.scale(2.6);
				c += b;
				//if( (*i).second == 1.0 ){
				glVertex3dv( &c[0] );
				//	}
				b.scale( -2.6);
				c += b;
				//	if( (*i).second == 1.0 ){
				glVertex3dv( &c[0] );
				// }

				glEnd();
			}

			// draw invis line then solid line
			for (int i = 0;i < 2;i++)
			{
				if (i == 0)
				{
					glDisable(GL_DEPTH_TEST);
					ot_set_rgb(ot_get_pref_rgb("loopcolor", OT3vec(1.0, 0.0, 1.0)), 0.6);
				}
				else
				{
					glEnable(GL_DEPTH_TEST);
					ot_set_rgb(ot_get_pref_rgb("loopcolor", OT3vec(1.0, 0.0, 1.0)), 1.0);
				}

				glBegin(GL_LINE_LOOP);
				for (vertlist_t::iterator i = b->traceverts.begin();
						i != b->traceverts.end(); ++i)
				{
					glVertex3d( (*i)->p.x, (*i)->p.y, (*i)->p.z );
				}
				glEnd();

				glPointSize(3.0);
				glBegin(GL_POINTS);
				for (vertlist_t::iterator i = b->traceverts.begin();
						i != b->traceverts.end(); ++i)
				{
					glVertex3d( (*i)->p.x, (*i)->p.y, (*i)->p.z );
				}
				glEnd();
			}
		}
	}

	glDisable(GL_TEXTURE_2D);
	// draw invis cursor then solid
	for (int i = 0;i < 2;i++)
	{
		if (i == 0)
		{
			glDisable(GL_DEPTH_TEST);
			glPointSize(3.0);
		}
		else
		{
			glEnable(GL_DEPTH_TEST);
			glPointSize(3.0);
		}
		ot_set_rgb(ot_get_pref_rgb("cursorcolor", OT3vec(0.0, 1.0, 0.0)), (i == 0) ? 0.4 : 1.0);

		glBegin(GL_POINTS);
		if (g_cursor_vis)
		{
			glVertex3d(g_cursor[0], g_cursor[1], g_cursor[2]);
		}
		glEnd();
	}

	glDisable(GL_DEPTH_TEST);

	/*
	glColor4f(0.0, 1.0, 0.0, 1.00);

	glBegin(GL_POINTS);
	int cc=0;
	for(vertlist_t::iterator i=g_traceverts.begin();
		i != g_traceverts.end(); ++i) {
			if(cc++==(countit%g_traceverts.size() )) glVertex3d( (*i)->p.x, (*i)->p.y, (*i)->p.z );
		}
	cc=0;
	glEnd();
	*/


	/*
	glColor4f(0.0, 0.0, 1.0, 1.00);
	glBegin(GL_POINTS);
	for(seglist_t::iterator i=tracesegs.begin(); i!=tracesegs.end(); ++i){
		if(cc==(countit )) glVertex3d( (*i)->v1->p.x, (*i)->v1->p.y, (*i)->v1->p.z );
		if(cc++==(countit )) glVertex3d( (*i)->v2->p.x, (*i)->v2->p.y, (*i)->v2->p.z );
	}
	glEnd();
	*/

	// draw 2d stroke
	ot_set_rgb(ot_get_pref_rgb("strokecolor", OT3vec(1.0, 0.0, 0.0)), 1.0);

	glPushMatrix();
	glLoadIdentity();

	glDisable(GL_DEPTH_TEST);
	glLineWidth(1.0f);
	glBegin(GL_LINES);

	int donecnt = g_points.size() - 1;
	for (int i = 0;i < donecnt;i++)
	{
		glVertex2d(g_points[i][0], g_points[i][1]);
		glVertex2d(g_points.wrap(i + 1)[0], g_points.wrap(i + 1)[1]);

	}
	glEnd();
	glBegin(GL_LINES);
	donecnt = g_over_points.size() - 1;
	for (int i = 0;i < donecnt;i++)
	{
		glVertex2d(g_over_points[i][0], g_over_points[i][1]);
		glVertex2d(g_over_points.wrap(i + 1)[0], g_over_points.wrap(i + 1)[1]);

	}
	glEnd();

	if (g_points.size() == 1)
	{
		glPointSize(3.0);

		glBegin(GL_POINTS);
		glVertex2d(g_points[0][0], g_points[0][1]);
		glEnd();
	}

	glColor3f(1.0, 1.0, 0.0);
	glPointSize(1.0);

	glBegin(GL_POINTS);
	//			glVertex2d(g_points[0][0], g_points[0][1]);
	//	glEnd();
	//glBegin(GL_LINES);
	donecnt = g_debugloop.size() - 1;
	for (int i = 0;i < donecnt;i += 2)
	{
		glVertex2d(g_debugloop[i][0], g_debugloop[i][1]);
		glVertex2d(g_debugloop.wrap(i + 1)[0], g_debugloop.wrap(i + 1)[1]);

	}
	glEnd();

	glColor3f(1.0, 1.0, 1.0);

	if (g_colorburst)
	{
		g_colorsample = true;
		// drop colorburst
		glBegin(GL_QUADS);
		glVertex2d( -1.0, -1.0);
		glVertex2d(1.0, -1.0);
		glVertex2d(1.0, 1.0);
		glVertex2d( -1.0, 1.0);
		glEnd();
	}

	if (g_colorsample)
	{
		glFlush();

		GLfloat readsample[3];
		glReadPixels(x, y, 1, 1, GL_RGB, GL_FLOAT, readsample);
		g_curcolor = OT3vec(readsample[0], readsample[1], readsample[2]);

		g_colorsample = false;
	}

	glEnable(GL_TEXTURE_2D);
	ot_display_status();

	glPopMatrix();


	glDepthMask(GL_TRUE);


	glFlush();
}


double ot_cursor_depth()
{
	if (!g_cursor_vis) return 0.5; // default
	OT3vec v;
	mproj(g_cursor[0], g_cursor[1], g_cursor[2], &v[0], &v[1], &v[2]);
	return v[2]; // depth
}

// extrude given loop starting at cursor z and moving up
GtsSurface* ot_extrude_loop_to_cursor(OT2loop *xloop, double d)
{
	// for now, just support boxes
	if (xloop->size() != 4) return NULL;

	double cursd = ot_cursor_depth();
	d = d * ot_measure_screen_at_depth(cursd);  // find global d value
	d /= ot_depth_range();  // convert do depth value

	GtsSurface* gtss = gts_surface_new (gts_surface_class (),
			GTS_FACE_CLASS (gts_face_class ()),
			GTS_EDGE_CLASS (gts_edge_class ()),
			GTS_VERTEX_CLASS (ot_gts_vertex_class ()));

	std::vector<GtsVertex*> frontvs, backvs;

	GtsVertex* v;
	// create front and back points
	for (int cycle = 0;cycle < 2;cycle++)
	{
		for (OT2loop::iterator i = xloop->begin(); i != xloop->end(); i++)
		{
			v = gts_vertex_new( gtss->vertex_class, 0.0, 0.0, 0.0);
			munproj_input( (*i)[0], (*i)[1],
					(cycle == 0) ? (cursd - d) : cursd, 	// cycle==0 is front, otherwise back
					&v->p.x, &v->p.y, &v->p.z);
			if (cycle == 0)
				frontvs.push_back(v);
			else
				backvs.push_back(v);
		}
	}
	// create side faces
	for (size_t i = 0;i < frontvs.size();i++)
	{
		GtsVertex* v[4];
		v[0] = frontvs[i]; v[1] = backvs[i];
		v[2] = frontvs[(i + 1) % frontvs.size()]; v[3] = backvs[(i + 1) % frontvs.size()];

		ot_gts_add_quad_from_existing_verts(gtss, v);
	}
	// create front faces
	{
		GtsVertex* v[4];
		v[0] = frontvs[0]; v[1] = frontvs[1]; v[2] = frontvs[3]; v[3] = frontvs[2];
		ot_gts_add_quad_from_existing_verts(gtss, v);

		v[1] = backvs[0]; v[0] = backvs[1]; v[3] = backvs[3]; v[2] = backvs[2];
		ot_gts_add_quad_from_existing_verts(gtss, v);
	}

	// face right way
	if (polygon_area(xloop) < 0.0 )
	{
		flip_faces(gtss);
	}
	return gtss;
}


void ot_main_reshape(GLint w, GLint h)
{
	// do standard reshaping
	win_w = w;
	win_h = h;
	glViewport(0, 0, w, h);

	setup_viewing_matrices();
}

static bool g_done = false;
int lastx, lasty;

// precision mode
bool g_pmode = false;

static gboolean ot_gts_refine_stop_number(gdouble cost, guint number, guint *max)
{
	if (number > *max)
		return TRUE;
	else return FALSE;
}

static long int g_framecnt;
static unsigned long g_tickcnt;

void ot_buttonisdown(double x, double y)
{
	static double lx = -1.0;
	static double ly = -1.0;

	if (ot_2d_dist(lx, ly, x, y) > 0.01)
	{
		lx = x;
		ly = y;
		addpt(x, y, &g_points, (g_over_points.size() > 0));
	}
}

bool g_dragging = false;
OTblobvec bsel;

void GLFWCALL OnMousePosChg(int x, int y)
{
	if (g_dragging)
	{
		/* drag in progress, simulate trackball */
		trackball(g_spin_quat,
				(2.0*lastx) / win_w - 1.0,
				( -2.0*lasty) / win_h + 1.0,
				(2.0*x) / win_w - 1.0,
				( -2.0*y) / win_h + 1.0
				 );
	}
	else
	{
		bool shiftheld = glfwGetKey(GLFW_KEY_LSHIFT) || glfwGetKey(GLFW_KEY_RSHIFT);
		if (glfwGetMouseButton(GLFW_MOUSE_BUTTON_LEFT) && !shiftheld)
		{
			double dx, dy;
			getmouse(0, 0, &dx, &dy);
			ot_buttonisdown(dx, dy);
		}
	}
}

void OnMouseBtnUp(int nBtn)
{
	if (g_dragging)
	{
		g_dragging = false;
		ot_view_setup_matrices(false);

		// if not a global transform, apply transform to selected parts
		if ( bsel.size() != g_allblobs.size() && (bsel.size() != 0) )
		{
			GLfloat m[4][4];
			build_rotmatrix(m, g_spin_quat);
			glLoadIdentity();
			glMultMatrixd(g_global_transform_inv);
			glMultMatrixf(&m[0][0]);
			glMultMatrixd(g_global_transform);

			double mat[16];
			glGetDoublev(GL_MODELVIEW_MATRIX, mat);

			for (OTblobvec::iterator i = bsel.begin(); i != bsel.end(); i++)
			{
				(*i)->matmult(mat);
			}
			ot_clear_all_modeltrans();

			// transform cursor as well
			glGetDoublev(GL_MODELVIEW_MATRIX, mat);
			g_cursor.matmult(mat);
		}
		else
		{
			GLfloat m[4][4];
			build_rotmatrix(m, g_spin_quat);
			ot_add_to_global_transformf(&m[0][0]);
		}
		clear_g_spin_quat();
	}
}

void OnMouseBtnDown(int nBtn)
{
	bool shiftheld = glfwGetKey(GLFW_KEY_LSHIFT) || glfwGetKey(GLFW_KEY_RSHIFT);

	g_debugloop.clear();

	double dx, dy;
	int x, y;
	getmouse(&x, &y, &dx, &dy);

	eventx = x;
	eventy = y;

	ot_clear_display_string();

	lastx = x;
	lasty = y;

	if (nBtn == GLFW_MOUSE_BUTTON_MIDDLE)
	{
		g_dragging = true;
		return ;
	}

	if (nBtn == GLFW_MOUSE_BUTTON_LEFT && !shiftheld && !g_dragging)
	{
		ot_buttonisdown(dx, dy);
	}

	if (nBtn == GLFW_MOUSE_BUTTON_RIGHT )
	{
		printf("gpoints size %d\n", g_points.size());
		if (g_points.size() < 1)
		{
			if (g_newselsurf != NULL)
			{
				// if shift
				//g_newselsurf->bSelected = !g_newselsurf->bSelected;
				// else
				g_newselsurf->bSelected = true;
			}
			g_cursor_vis = true;
			g_cursor = g_newcursor;

			return ;
		}


		if (bsel.size() <= 0)
		{
			if (ot_polygon_intersecting(&g_points))
			{
				ot_set_display_string("Stroke is self-intersecting.");
				return ;
			}

			GtsSurface* newsurf = NULL;
			if (g_points.size() < 5)
			{
				// do special things for low number of points
				OT2loop xloop;
				OT2vec depthvec;
				if (g_points.size() == 3)
				{
					// topleft, bottom right, then depth
					xloop.push_back(g_points[0]);
					xloop.push_back(OT2vec(g_points[0][0], g_points[1][1]));
					xloop.push_back(g_points[1]);
					xloop.push_back(OT2vec(g_points[1][0], g_points[0][1]) );
					depthvec = g_points[2];
					depthvec -= g_points[1];
				}
				if (g_points.size() == 4)
				{
					// topleft, bottom right, then depth
					xloop.push_back(g_points[0]);
					xloop.push_back(OT2vec(g_points[0][0], g_points[2][1]));
					xloop.push_back(g_points[2]);
					xloop.push_back(OT2vec(g_points[2][0], g_points[0][1]) );
					depthvec = g_points[3];
					depthvec -= g_points[2];
				}

				if (xloop.size()) newsurf = ot_extrude_loop_to_cursor(&xloop, depthvec.len());
				if (!newsurf)
				{
					ot_set_display_string("Not enough vertices to extrude.");
					return ;
				}
			}

			if (!newsurf)
			{
				if (!g_pmode)
					ot_loop_simplify(&g_points, true);

				// convert format for input to inflater
				double* d = new double[g_points.size() * 2];
				for (size_t i = 0;i < g_points.size();i++)
				{
					d[i*2 + 0] = g_points[i][0];
					d[i*2 + 1] = g_points[i][1];
				}

				g_opi = new ot_public_inflate (d, g_points.size(), ot_get_pref_int("edgetess", 4) );

				delete[] d;

				newsurf = gts_surface_new (gts_surface_class (),
						GTS_FACE_CLASS (gts_face_class ()),
						GTS_EDGE_CLASS (gts_edge_class ()),
						GTS_VERTEX_CLASS (ot_gts_vertex_class ()));
				newsurf = g_opi->generate_gts(newsurf);

				// change to current view
				ot_gts_vertices_to_global(newsurf, ot_cursor_depth() /* desired depth */ );
				ot_calculate_geometry_warning(newsurf);

				// smooth it
				gdouble maxcosine2 = cos ( ( /*deg*/90.0 * M_PI) / 180.);
				maxcosine2 *= maxcosine2;

// 				if (ot_get_pref_int("auto_laplace", 1))
// 				{
// 					ex_gts_laplace_filter(newsurf, ot_get_pref_int("laplace_niter", 5),
// 							ot_get_pref_double("laplace_lambda", 0.1),
// 							maxcosine2,  /*fold*/false);
// 				}

				double auto_simplify = ot_get_pref_double("auto_simplify", 0.5);
				if (auto_simplify > 0.0)
				{
					guint number = (gts_surface_edge_number(newsurf) * auto_simplify);
					/* "dumb" coarsen		*/
					gts_surface_coarsen(newsurf, NULL,
							NULL, NULL, NULL,
							(GtsStopFunc) gts_coarsen_stop_number, &number,
							M_PI / 180.0 //minangle
									   );
				}
			}

			OTblob* b = new OTblob();
			b->gtss = newsurf;
			ot_calculate_geometry_warning(b->gtss);

			g_points.clear();

			g_allblobs.push_back(b);
		}
		else
		{
			if ( bsel.size() != 1)
			{
				ot_set_display_string("Only one surface may be selected.");
				return ;
			}
			OTblob* b = (* bsel.begin() );

			if (b->traceverts.size() > 0)
			{
				if (!g_pmode)
					ot_loop_simplify(&g_points, false, 2.0);

				perform_action(OT_PERFORM_PUNCHOUT, b);
				extrude_trace(&g_points, &b->traceverts, b);
				b->traceverts.clear();

				g_points.clear();
			}
			else
			{
				if (ot_polygon_intersecting(&g_points))
				{
					ot_set_display_string("Stroke is self-intersecting.");

					return ;
				}

				if (!g_pmode)
					ot_loop_simplify(&g_points, true, 2.0);
				ot_realize_loop_on_surface(&g_points, &b->traceverts, b, false);

				g_points.clear();
			}
		}
	}
	else
	{
		/*		static REAL lx=-1.0;
				static REAL ly=-1.0;
				REAL x,y;
				getmouse(x,y);
				if(dist(lx, ly, x, y) > 0.01){
					lx = x;
					ly = y;
					addpt(x,y);
				}
		*/
	}
}

void GLFWCALL OnMouseBtnChg(int nBtn, int nType)
{
	if (nType == GLFW_PRESS)
	{
		OnMouseBtnDown(nBtn);
	}
	else if (nType == GLFW_RELEASE)
	{
		OnMouseBtnUp(nBtn);
	}
}

static int do_coarsen = 0;

static bool repeat_last = false;
static int pklast = 0;

void GLFWCALL OnKeyChg(int nkey, int action)
{
	ot_selblob(&g_allblobs, &bsel);
	if (action != GLFW_PRESS )
	{
		return ;
	}

	bool shiftdown = glfwGetKey(GLFW_KEY_LSHIFT) || glfwGetKey(GLFW_KEY_RSHIFT);
	bool controldown = glfwGetKey(GLFW_KEY_LCTRL) || glfwGetKey(GLFW_KEY_RCTRL);

	if (nkey == GLFW_KEY_F2) draw_2d = !draw_2d;
	// TEMP counter,  and place cursor
	if (nkey == GLFW_KEY_SPACE)
	{

		g_cursor_vis = true;
		if (shiftdown)
			g_cursor = g_defaultcursor;
		else
			g_cursor = g_newcursor;

		countit++;
	}

	// laplace
// 	if ( nkey == GLFW_KEY_SPACE )
// 	{
// 		// smooth it
// 		gdouble maxcosine2 = cos ( ( /*deg*/90.0 * M_PI) / 180.);
// 		maxcosine2 *= maxcosine2;
// 
// 		for (OTblobvec::iterator i = bsel.begin(); i != bsel.end(); i++)
// 		{
// 			ex_gts_laplace_filter((*i)->gtss, ot_get_pref_int("laplace_manual_niter", 1),
// 					ot_get_pref_double("laplace_manual_lambda", 0.1),
// 					maxcosine2,  /*fold*/false);
// 		}
// 	}

	// cut out
	if ( nkey == 'X' )
	{
		for (OTblobvec::iterator i = bsel.begin(); i != bsel.end(); i++)
		{
			if (ot_realize_loop_on_surface(&g_points, &(*i)->traceverts, (*i), true))
			{
				if (controldown)
				{
					perform_action(OT_PERFORM_PERFORATE, (*i));
				}
				else if (shiftdown)
				{
					perform_action(OT_PERFORM_INV_CUT, (*i));
				}
				else
				{
					perform_action(OT_PERFORM_CUT, (*i));
				}

				g_points.clear();
			}
			else
			{
				ot_set_display_string("Stroke does not intersect surface.");
			}
		}
	}

	// select all
	if ( nkey == 'A' )
	{
		bool anysel = false;
		for (OTblobvec::iterator i = g_allblobs.begin(); i != g_allblobs.end(); i++)
		{
			(*i)->bSelected = true;
		}

		return ;
	}
	// unselect
	if ( nkey == GLFW_KEY_ESC )
	{
		if (g_points.size())
		{
			g_points.clear();
			return ;
		}

		bool anysel = false;
		for (OTblobvec::iterator i = g_allblobs.begin(); i != g_allblobs.end(); i++)
		{
			OT_ASSERT( (*i) );
			anysel = anysel || (*i)->bSelected;
			(*i)->bSelected = false;
		}
		if (anysel) return ;
		g_cursor_vis = false;

		return ;
	}
	if ( nkey == GLFW_KEY_UP )
	{
		g_scale *= 1.5;
	}
	if ( nkey == GLFW_KEY_DOWN )
	{
		g_scale *= 1.0 / 1.5 ;
		if (g_scale == 0.0) g_scale = 1e-6;
	}
	if (g_scale != 1.0)
	{
		g_scale_actual *= g_scale;

		glMatrixMode(GL_MODELVIEW);
		glLoadIdentity();
		glScalef(1.0 / g_scale, 1.0 / g_scale, 1.0 / g_scale);
		glMultMatrixd(g_global_transform);

		ot_modelview_to_global_transform();
		g_scale = 1.0;
	}

	// simplify
	if ( nkey == 'I' )
	{
		do_coarsen = shiftdown ? 0 : 1;
		if (controldown) do_coarsen += 2;

		ot_set_display_string("Simplifying selected surfaces...");
		do_simplify();
	}

	// delete selected or if animation, delete point
	if ( nkey == GLFW_KEY_BACKSPACE )
	{
		for (OTblobvec::iterator i = g_allblobs.begin(); i != g_allblobs.end();)
		{
			OT_ASSERT(*i);
			bool nexti = true;
			if ( (*i)->bSelected )
			{
				delete (*i);
				i = g_allblobs.erase(i);
				nexti = false;
			}
			if (nexti) i++;
		}
	}

	// hide
	if ( nkey == 'H' )
	{
		if (shiftdown)
			for (OTblobvec::iterator i = g_allblobs.begin(); i != g_allblobs.end(); i++)
			{
				(*i)->invis = false;
			}
		else
		{
			for (OTblobvec::iterator i = bsel.begin(); i != bsel.end(); i++)
			{
				(*i)->invis = true;
				(*i)->bSelected = false;
			}
		}
	}

	// new origin
	if ( nkey == 'O' )
	{
		if (g_cursor_vis)
		{
			glMatrixMode(GL_MODELVIEW);
			glLoadIdentity();
			OT3vec oldcursor = g_cursor;
			g_cursor.matmult(g_global_transform);
			glTranslatef( -g_cursor[0], -g_cursor[1], -g_cursor[2]);
			g_cursor = oldcursor;
			glMultMatrixd(g_global_transform);

			ot_modelview_to_global_transform();
		}
		if (shiftdown)
		{
			/* permanent movement */
			for (OTblobvec::iterator i = g_allblobs.begin(); i != g_allblobs.end(); i++)
			{
				ot_gts_new_origin( (*i) );
			}

			// transform cursor as well
			double mat[16];
			glGetDoublev(GL_MODELVIEW_MATRIX, mat);
			g_cursor.matmult(mat);

			// reset global transform
			glMatrixMode(GL_MODELVIEW);
			glLoadIdentity();
			ot_modelview_to_global_transform();
		}
	}

	//join
	if ( nkey == 'J' )
	{
		if (bsel.size() <= 1)
		{
			ot_set_display_string("Join requires more than one surface selected.");
			return ;
		}
		OTblob* b = NULL;
		for (OTblobvec::iterator i = g_allblobs.begin(); i != g_allblobs.end();)
		{
			OT_ASSERT(*i);
			bool nexti = true;
			if ( (*i)->bSelected )
			{
				if (b == NULL)
				{
					b = *i;
				}
				else
				{
					if (b->join( *i ))
					{
						delete (*i);
						i = g_allblobs.erase(i);
						nexti = false;
					}
				}
			}
			if (nexti) i++;
		}

	}
	//precision mode
	if ( nkey == 'P' )
	{
		g_pmode = !g_pmode;
		if (g_pmode)
		{
			ot_set_display_string("Precision mode on.");
		}
		else
		{
			ot_set_display_string("Precision mode off.");
		}
	}

	if ( nkey == GLFW_KEY_TAB )
	{
		g_auto_select = !g_auto_select;
		if (g_auto_select)
		{
			ot_set_display_string("Auto-select mode on.");
		}
		else
		{
			ot_set_display_string("Auto-select mode off.");
		}
	}



	// select
	if ( nkey == 'E' )
	{
		if (shiftdown){ g_pickface = 1; return ;}
		if (controldown){ g_pickface = 2; return ;}

		for (OTblobvec::iterator i = bsel.begin(); i != bsel.end(); i++)
		{
			OT_ASSERT(*i);
			(*i)->selsel();
		}
	}
	// oversketch
// 	if ( nkey == 'V' )
// 	{
// 		if (g_points.size() < 2) return ;
// 		if (g_over_points.size() == 0)
// 		{
// 			// this is here for conveience- if nothing is selected, assume we mean to do it on everything
// 			for (OTblobvec::iterator i = bsel.begin(); i != bsel.end(); i++)
// 			{
// 				OT_ASSERT(*i);
// 				if ( (*i)->vertsel.size() == 0)
// 					(*i)->selsel();
// 			}
// 
// 			g_over_points = OT2loop(g_points);
// 			g_points.clear();
// 			return ;
// 		}
// 		ot_loop_regulate(&g_points, false, 22);
// 		ot_loop_regulate(&g_over_points, false, 22);
// 
// 		for (OTblobvec::iterator i = bsel.begin(); i != bsel.end(); i++)
// 		{
// 			OT_ASSERT(*i);
// 			(*i)->oversketch(&g_over_points, &g_points);
// 		}
// 
// 		g_points.clear(); g_over_points.clear();
// 	}
	// simplify loop
	if ( nkey == 'Z' )
	{
		ot_loop_simplify(&g_points, false, 2.0);
	}

	// render mode
	if ( nkey == 'r' )
	{
		if (shiftdown)
		{
			g_show_wireframe = !g_show_wireframe;
		}
		else
		{
			g_render_mode = (render_mode_t)( ( (int) g_render_mode + 1) % ((int)OT_RENDER_SOLID + 1) );
			if (g_render_mode == OT_RENDER_WIREFRAME)
				g_show_wireframe = true;
			else
				g_show_wireframe = false;
		}
	}

	// select color
	if ( nkey == 'C' )
	{
		if (shiftdown){ g_colorsample = true; }
		else
			g_colorburst = !g_colorburst;
		return ;
	}

}

bool do_simplify()
{
	ot_selblob(&g_allblobs, &bsel);

	if (bsel.size() <= 0)
	{
		ot_set_display_string("Simplify requires at least one surface selected.");
		return false;
	}

	int tv = 0, tf = 0, te = 0;
	for (OTblobvec::iterator i = bsel.begin(); i != bsel.end(); i++)
	{
		OT_ASSERT(*i);
		GtsSurface* gtss = (*i)->gtss;
		(*i)->clear_si();
		(*i)->traceverts.clear();

		/* WARNING, CAN ONLY DO THIS WHEN THERE IS NO TRACE ON SURFACE,
			OR ELSE TRACE POINTS MAY BECOME INVALID 
		GtsVolumeOptimizedParams voparams = { 0.5, 0.5, 0.0 };
		gdouble cmax = 1e-12;

		if( (traceverts->size() == 0) && gtss) gts_surface_coarsen(gtss, (GtsKeyFunc) gts_volume_optimized_cost,
			&voparams, (GtsCoarsenFunc) gts_volume_optimized_vertex, &voparams,
			(GtsStopFunc) gts_coarsen_stop_cost, &cmax,
			0.01 //minangle
			);*/ 
		/* alternate simp by count */
		// eliminate useless triangles
		GtsVolumeOptimizedParams voparams = { 0.5, 0.5, 0. };

		if (do_coarsen == 1)
		{
			guint number = (gts_surface_edge_number(gtss) * 9) / 10;
			if (number > 6)
			{
				//					number = ot_get_pref_int("edgemax", 1000);

				gts_surface_coarsen(gtss, (GtsKeyFunc) gts_volume_optimized_cost,
						&voparams, (GtsCoarsenFunc) gts_volume_optimized_vertex, &voparams,
						(GtsStopFunc) gts_coarsen_stop_number, &number,
						M_PI / 180.0 //minangle
								   );
			}
		}
		else if (do_coarsen == 0)
		{
			guint number = (gts_surface_edge_number(gtss) * 130) / 100;
			gts_surface_refine (gtss,
					NULL /* cost func*/,
					NULL /* params */,
					NULL /* refine func*/,
					NULL /* refine data */,
					(GtsStopFunc) ot_gts_refine_stop_number, &number);
		}
		else
		{
			guint number = (gts_surface_edge_number(gtss) * 9) / 10;
			/* "dumb" coarsen		*/
			gts_surface_coarsen(gtss, NULL,
					&voparams, NULL, &voparams,
					(GtsStopFunc) gts_coarsen_stop_number, &number,
					M_PI / 180.0 //minangle
							   );
		}

		te += gts_surface_edge_number(gtss);
		tv += gts_surface_vertex_number(gtss);
		tf += gts_surface_face_number(gtss);
	}
	static char displaystr[256];
	sprintf(displaystr, "Done.  Faces:%d Edges:%d Verts:%d", tf, te, tv);
	ot_set_display_string(displaystr);
}

int main(int argc, char *argv[])
{
	ot_read_prefs();
	g_auto_select = true;

	// Initialise GLFW
	glfwInit();

	// Open OpenGL window
	if ( !glfwOpenWindow(win_w = 800, win_h = 600, 0, 0, 0, 0, 0, 0, GLFW_WINDOW ) )
	{
		glfwTerminate();
		return 0;
	}

	// Set tile
	glfwSetWindowTitle( "OpenTeddy" );

	// Set callback function
	glfwSetKeyCallback( OnKeyChg );
	glfwSetMousePosCallback(OnMousePosChg);
	glfwSetMouseButtonCallback(OnMouseBtnChg);
	glfwSetWindowSizeCallback(ot_main_reshape);

	ot_main_reshape(win_w, win_h);
	printf("GL_RENDERER   = %s\n", (char *)	glGetString(GL_RENDERER));
	printf("GL_VERSION    = %s\n", (char *) glGetString(GL_VERSION));
	printf("GL_VENDOR     = %s\n", (char *) glGetString(GL_VENDOR));
	printf("GL_EXTENSIONS = %s\n", (char *)	glGetString(GL_EXTENSIONS));

	// Main loop
	g_done = GL_TRUE;
	while ( g_done )
	{
		tickdraw();

		// Swap buffers
		glfwSwapBuffers();

		// Check if the window was closed
		g_done = g_done && glfwGetWindowParam( GLFW_OPENED );
		glfwSleep(0.005);
	}

	// Close OpenGL window and terminate GLFW
	glfwTerminate();

	return 0;
}